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ABSTRACT 


Emitter  location  finding  enables  important  functionality  for  both  military  and 
civilian  applications.  GPS  is  the  most  recognized  and  widely  used  positioning 
system,  but  it  is  a  receiver  location  system  that  functions  in  a  markedly  different 
manner  from  emitter  location.  Many  geo-location  techniques  predate  and  have 
been  proposed  as  an  alternative  to  GPS.  Some  of  the  more  commonly  used  and 
exploited  of  these  techniques  are  angle  of  arrival,  time  of  arrival,  frequency 
difference  of  arrival,  and  time  difference  of  arrival  (TDOA).  This  thesis  is  primarily 
focused  on  TDOA. 

These  techniques  are  important  for  military  applications.  Location  finding 
is  a  part  of  electronic  warfare  support,  which  is  one  of  the  main  braches  of 
electronic  warfare.  Because  these  techniques  are  platform  independent,  they  can 
be  used  with  any  system  or  platform,  such  as  UAVs,  manned  aircraft,  ground 
locations,  etc.  In  Turkey  it  is  vitally  important  for  the  army  conducting  search  and 
destroy  operations  against  terrorists  to  locate  emitters  associated  with  these 
terrorists. 

The  simulation  developed  in  this  thesis  provides  a  better  understanding  of 
the  accuracy  of  TDOA  based  geolocation  systems.  Combinations  of  receivers 
and  techniques  are  explored  to  determine  the  optimal  solutions.  The  factors  of 
noise  and  distance  have  a  linear  effect  on  accuracy.  The  best  combination  of 
receivers  is  determined  with  consideration  to  using  a  combination  of  fixed  and 
airborne  platforms.  The  best  distribution  for  highest  accuracy  is  determined  and 
discussed. 
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I.  INTRODUCTION 


A.  AREA  OF  RESEARCH 

Electronic  warfare  (EW)  plays  a  dominant  role  in  today’s  world  of 
technological  warfare.  It  has  been  considered  to  be  a  force  multiplier  for 
decades.  Nations  who  understand  the  importance  of  the  EW  have  always  used  it 
and  benefited.  We  as  the  soldiers  cannot  imagine  a  battlefield  without  EW  and  its 
components.  “EW  is  one  of  the  five  core  capabilities”  (JP  3-13-1,  2007).  Joint 
Publication  (JP)  3-13-1  Electronic  Warfare  is  U.S.  joint  doctrine  for  EW  and 
provides  joint  doctrine  for  every  part  of  EW,  ranging  from  planning  and 
preparation  to  execution  for  military  operations.  This  publication  is  an  important 
source  for  anyone  wanting  to  define  what  EW  is  and  what  components  compose 
it.  We  cannot  win  a  war  or  even  a  battle  without  the  use  of  EW. 

EW  includes  three  major  branches:  Electronic  Attack  (EA),  Electronic 
Protection  (EP)  and  Electronic  Warfare  Support  (ES).  “EA  involves  the  use  of  EM 
energy,  directed  energy,  or  anti-radiation  weapons  to  attack  personnel,  facilities, 
or  equipment  with  the  intent  of  degrading,  neutralizing,  or  destroying  enemy 
combat  capability  and  is  considered  a  form  of  fire.  EP  involves  actions  taken  to 
protect  personnel,  facilities,  and  equipment  from  any  effects  of  friendly  or  enemy 
use  of  the  electromagnetic  spectrum  that  degrade,  neutralize,  or  destroy  friendly 
combat  capability.  ES  is  the  subdivision  of  EW  involving  actions  tasked  by,  or 
under  direct  control  of,  an  operational  commander  to  search  for,  intercept, 
identify,  and  locate  or  localize  sources  of  intentional  and  unintentional  radiated 
EM  energy  for  the  purpose  of  immediate  threat  recognition,  targeting,  planning, 
and  conduct  of  future  operations”  (JP  3-1 3-1 , 2007). 

ES  includes  direction  finding  of  enemy  RF  signals.  Direction  finding  can 
also  be  considered  as  finding  the  enemy  location  by  analyzing  the  signal  or 
signals  that  are  transmitted  by  the  emitter.  This  is  also  called  emitter  location  or, 
if  the  emitter  is  on  the  earth’s  surface,  emitter  geo-location. 
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It  is  important  for  both  the  military  and  civilians  to  be  able  to  find  the 
location  of  RF  emitters  on  the  surface  of  the  earth  under  extreme  conditions, 
such  as  in  war  for  the  military  and  in  the  case  of  accidents  or  disasters  for  civilian 
response.  Military  forces  want  to  know  the  where  the  enemy  is.  Determining  the 
location  of  enemy  forces  has  always  been  important  to  combat,  and  various 
methods  have  been  studied  and  developed  through  the  ages.  The  easiest  and 
most  accurate  way  to  do  this  historically  has  been  to  send  a  man  forward  to 
locate  the  enemy.  Time  and  available  human  resources  limit  this.  There  are 
inherent  risks,  such  as  having  the  soldiers  killed  and  being  recognized  by  the 
enemy  forces  if  they  capture  the  reconnaissance  team.  Because  of  these 
limitations  people  have  tried  to  use  technology  as  much  as  possible  in  the 
modern  era. 

It  is  important  to  know  the  position  of  the  enemy  without  letting  them  know 
that  you  are  trying  to  locate  them.  This  can  be  done  with  passive  location  finding 
techniques.  In  contrast,  people  have  been  using  active  location  finding 
techniques,  such  as  radar,  for  decades.  When  active  location  finding  techniques 
are  used  they  let  the  enemy  know  that  they  are  being  identified  or  located 
because  active  location  finding  techniques  uses  detectable  RF  signals. 

To  reduce  the  chance  that  the  enemy  will  discover  that  they  are  under 
scrutiny,  scientists  and  engineers  have  developed  techniques  for  passive 
location  finding.  These  techniques  depend  upon  the  interception  of  enemy 
transmissions  within  the  EM  spectrum,  rather  than  emanating  friendly  EM 
transmissions.  Some  of  the  more  common  of  these  techniques  are 

•  TDOA  (Time  Difference  of  Arrival) 

•  FDOA  (Frequency  Difference  of  Arrival) 

•  AOA  (Angle  of  Arrival) 

•  TOA  (Time  of  Arrival) 
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Among  these  techniques,  TDOA  is  currently  the  most  commonly  studied 
technique.  This  technique  is  platform  independent.  TDOA  is  most  easily  applied 
to  pulsed  signals.  Since  all  signals  of  interest  are  not  pulsed  signals,  some 
techniques  have  to  be  implemented  for  continuous  wave  (CW)  signals.  Pulsed 
signals  have  distinctive  features  that  can  be  recognized  in  the  waveform.  CW 
signals  do  not  have  the  same  distinctive  features  and  can  change  over  time.  The 
signal  depends  on  the  message  that  is  carried  over  the  carrier.  Since  the 
message  over  the  carrier  is  random,  it  is  difficult  to  determine  the  beginning  and 
the  end  of  the  signal. 

Radars  typically  use  pulsed  signals,  and  communication  systems  typically 
use  various  modulations  of  CW  signals.  Ground  forces  are  normally  more 
interested  in  communication  signals  than  radar;  military  forces  need  to  coordinate 
their  needs  with  scientists  to  find  an  appropriate  geolocation  technique  for  the 
CW  signals  associated  with  these  communications. 

The  FDOA  technique  requires  a  long-duration  signal  to  be  able  to 
determine  direction  to  a  sufficient  accuracy  for  position  determination. 
Conversely,  the  TDOA  technique  does  not  need  this  long  duration.  Instead  of 
requiring  a  long  duration  signal,  very  short  duration  signals  can  be  used  to 
determine  the  location  if  all  the  necessary  receivers  receive  the  signal.  Proper 
coordination  of  the  signals  between  the  receivers  and  the  main  node  that 
compute  the  position,  and  synchronization  among  the  receivers,  are  necessary 
for  accurate  emitter  position  estimation. 

For  ground  forces,  passive  location  finding  becomes  important  for  locating 

enemy  combatants  or  terrorists  in  mountainous  areas,  such  as  those  found  in  the 

southeast  part  of  Turkey.  Turkish  land  forces  have  been  conducting  search  and 

destroy  operations  for  over  30  years  in  this  region.  As  a  part  of  these  operations, 

Turkish  Land  Forces  typically  first  attempt  to  locate  the  terrorists  who  hide  within 

the  cave  infrastructure  common  to  the  area  and  then  destroy  them.  Terrorists  in 

the  region  are  known  to  move  as  small  groups  and  to  use  frequency  modulation 

(FM)  based  communication  devices.  If  friendly  forces  can  use  a  passive 
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technique  to  locate  them  this  would  make  the  search  process  easier.  Depending 
upon  the  operation,  an  accuracy  of  100-200  meters  is  sufficient.  For  this 
purpose,  unmanned  air  vehicles  (UAVs)  can  be  used.  Since  typical  insurgent  FM 
communications  are  not  long-duration  or  long  distance,  troops  should  have 
mobile  direction  finding  (DF)  devices  with  them  or  nearby. 

High  clutter  brings  out  another  problem.  When  troops  use  DF  devices  in 
mountainous  areas  such  as  southeast  Turkey,  devices  might  receive  signals 
more  than  once  due  to  the  reflection  of  the  signal  from  the  clutter.  This  makes  the 
problem  more  complicated  and  more  difficult,  especially  when  troops  want  to 
locate  the  enemy  emitter  in  three  dimensions. 

1.  Tactical  Situation 

This  thesis  is  concerned  with  a  tactical  situation  in  Turkey. 

Turkish  armed  forces  has  been  planning  and  conducting  military 
operations  against  the  PKK  (Partiya  Karker  Kurdistan  (Kurdistan  Worker's 
Party)).  The  PKK,  whose  main  purpose  is  to  weaken  the  ruling  government  in 
Turkey,  is  classified  as  a  terrorist  organization.  Most  nations  of  the  world  accept 
this  classification. 

These  operations  are  conducted  mostly  in  the  southeast  region  of  the 
Turkey.  That  part  of  Turkey  is  very  mountainous  and  includes  many  caves. 
Terrorists  hide  in  the  caves  and  come  out  whenever  they  want  to  attack. 
Operating  within  mountainous  terrain  gives  them  an  advantage  in  hiding  from 
government  forces  and  helps  them  to  move  more  discretely.  If  military  forces  try 
to  follow  them,  especially  at  night,  they  escape  easily  because  they  know  the 
terrain  better  than  the  military  forces.  Darkness  covers  them  and  mountains 
hinder  the  movement  of  the  larger  military  forces  more  than  the  terrorists. 

These  terrorists  mostly  come  from  Iraq,  Iran  and  some  other  neighbor 
countries  where  they  receive  their  training.  They  learn  how  to  fight  and  use 
mines,  explosives  and  other  weapons  and  equipment.  They  cross  the  border  in 
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small  groups  at  night.  They  travel  to  their  destinations  on  foot,  sometimes  taking 
from  one  to  two  weeks. 

The  terrorists  typically  use  hand-held  radio  devices  for  communication 
between  groups  and  the  main  camp  if  they  are  close  to  it.  Some  high-ranking 
terrorists  use  cell  phones  or  satellite  phones  if  they  can,  but  GSM  companies  do 
not  cover  most  of  the  mountainous  areas  where  they  operate.  The  devices  used 
between  groups  normally  use  FM  modulation.  They  try  to  use  these  devices  as 
little  as  possible  to  avoid  detection  and  location.  Often  the  duration  of  these 
communications  number  in  the  seconds  or  at  best  tens  of  seconds. 

FDOA  cannot  be  used  because  of  its  long  duration  signal  requirement. 
TDOA  therefore,  since  it  does  not  have  the  same  requirement,  might  be  the  best 
solution. 

Turkish  forces  have  been  using  some  technology  and  intelligence  support 
to  find  the  terrorists,  but  they  still  manage  to  escape  using  the  methods 
discussed  above. 

These  terrorists  must  be  found  and  destroyed.  A  better  method  to  find 
them  is  to  use  some  sort  of  geo-location  system  targeted  against  their 
communications  signals.  Even  though  they  try  to  use  FM  as  little  as  possible, 
TDOA  is  a  very  good  technique  for  this  purpose.  If  land  forces  can  use  a  tactic 
based  on  TDOA  they  will  be  more  successful  in  finding,  locating,  and  destroying 
the  groups  of  terrorists. 

2.  TDOA 

Time  Difference  of  Arrival  is  one  of  the  basic  DF  techniques  used  to  locate 
RF  emitters  for  CW  signals.  “TDOA  takes  advantage  of  the  fact  that  a  transmitted 
signal  will  arrive  at  different  sensors  at  different  times.”  (Batson,  McEachen,  & 
Tummala,  2012)  “A  number  of  spatially  separated  sensors  capture  the  emitted 
signal  and  the  time  differences  of  arrival  (TDOAs)  at  the  sensors  are  determined. 


5 


Using  the  TDOA’s,  emitter  location  relative  to  the  sensor  can  be  calculated.” 
(Chan  &  Ho,  1994) 

B.  MAJOR  RESEARCH  QUESTIONS 

This  thesis  will  explain  the  following  subjects. 

•  What  are  EW  and  its  components? 

•  What  is  the  role  of  DF  in  EW? 

This  study  will  answer  the  following  questions: 

•  What  are  the  advantages  and  disadvantages  of  using  TDOA? 

•  How  effectively  can  ground  troops  use  TDOA  DF  techniques  in  a 
high  clutter  area  with  noise? 

•  How  can  we  simulate  TDOA  technique? 

•  What  are  the  optimal  uses  of  the  TDOA  technique? 

C.  IMPORTANCE  AND  BENEFITS  OF  THE  STUDY 

This  study  can  be  used  as  a  guide  for  developing  a  geolocation  system 
based  on  TDOA.  It  will  also  identify  the  benefits  of  using  a  TDOA-based 
geolocation  system  for  ground  forces  in  high  clutter  terrain.  It  will  clarify  the  major 
EW  components  and  their  usage  as  an  introduction  to  discussions  of  TDOA. 

D.  ORGANIZATION  OF  THE  THESIS 

This  thesis  is  composed  of  five  chapters.  Chapter  I  introduces  the  area  of 
research,  major  research  questions,  importance  and  benefits  of  the  study  and  the 
organization  of  the  thesis. 

Chapter  II  presents  information  about  EW  fundamentals.  EW  and  its 
components  are  examined.  The  relationship  between  EW  and  10  is  explained. 
The  importance  of  EW  is  explained. 
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Chapter  III  explains  DF  techniques  such  as  directional  antenna,  phase 
interferometry,  TDOA,  and  FDOA.  TDOA,  as  a  part  of  ES,  is  explained  briefly. 
The  details  of  the  TDOA  technique  and  its  usage  in  the  field  are  examined. 
Capabilities  of  TDOA  are  described.  The  developed  TDOA-based  location 
system  for  ground  forces  in  a  high  clutter  area  are  explained  from  Ezzat’s  paper. 
The  least  square  estimation  is  explained  and  the  required  recursive  least  square 
estimation  approach  is  examined  in  detail.  Sources  of  error  are  explained  for  DF 
and  TDOA. 

Chapter  IV  explains  the  developed  simulation  for  the  system.  Assumptions 
and  restrictions  for  the  simulation  are  listed.  After  running  the  simulation,  five 
types  of  analyses  are  done  to  see  their  effects  on  the  accuracy  of  the  proposed 
system  and  the  simulation. 

Chapter  V  is  the  conclusion.  The  results  obtained  from  the  system 
simulation  are  explained.  Recommendations  for  future  works  are  given. 
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II.  ELECTRONIC  WARFARE 


In  this  chapter,  information  operations  and  its  relationship  to  electronic 
warfare  are  defined.  Then  electronic  warfare  and  its  components  explained. 
Since  the  topic  of  the  thesis  is  direction  finding,  and  because  direction  finding  is  a 
subdivision  of  electronic  warfare  support  (ES),  ES  is  explained  in  more  detail. 
These  definitions  are  made  according  to  Frater  and  Ryan’s  book  on  EW,  with 
secondary  reference  using  U.S.  EW  doctrine  as  discussed  in  Joint  Publication 
3.13.1. 

Another  classification  of  EW  components  is  based  on  their  activities, 
which  can  be  defined  as  either  active  or  passive.  Active/passive  classification  is 
discussed. 

A.  INFORMATION  OPERATIONS 

Information  has  a  vital  importance  for  any  nation  and  for  its  existence. 
Every  nation  around  the  world  tries  to  obtain  information,  use  it  as  much  as 
possible  for  its  own  good,  and  prevent  the  enemy  from  taking  advantage  of  it. 
Information  operations  (10)  is  the  effort  done  for  this  purpose,  by  human  or  by 
machine.  10  incorporates  the  actions  taken  to  preserve  the  integrity  of  one’s  own 
information  system  from  exploitation,  corruption,  or  disruption,  while  at  the  same 
time  exploiting,  corrupting  and  destroying  the  adversary’s  information  system 
(Adamy,  2004). 

Using  U.S.  doctrine  as  a  reference,  10  coordinates  and  synchronizes  five 
core  capabilities  to  help  the  commander  to  reach  his/her  purpose  in  the 
battlefield.  These  capabilities  are  psychological  operations,  military  deception, 
operations  security,  electronic  warfare  and  computer  network  operations  (JP 
3.13,  2006).  In  addition  to  the  core  capabilities  doctrine  defines  supporting 
capabilities,  which  are  information  assurance  (IA),  physical  security,  physical 
attack,  counterintelligence,  and  combat  camera,  (JP  3-13,  2006)  and  related 
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capabilities  of  civil  military  operations,  public  affairs  and  defense  support  to  public 
diplomacy  (JP  3-13,  2006). 

The  world  has  gone  to  cellular  and  wireless.  Since  the  Information  Age 
produced  a  revolution  in  military  operations,  the  electromagnetic  spectrum  is 
essential  to  the  transmission  and  reception  of  information  in  the  modern  world, 
and  therefore  EW  is  a  critical  element  of  information  operations.  As 
communication  and  information  systems  become  increasingly  vital  for  military 
and  civilian  society,  they  can  become  targets  in  war  for  an  enemy;  therefore, 
they  can  play  a  significant  role  for  offensive  and  defensive  operations.  The 
military  has  adopted  communication  and  information  systems  and  they  have 
become  essential  for  military  operations.  Commanders  need  information  to  be 
able  to  cope  with  the  complexity  of  modern  warfare.  This  reliance  on  information 
in  turn  makes  them  vulnerable  to  attack.  There  is  an  emphasis  for  commanders 
to  attack  adversary  information  and  communication  systems.  Modern  battlefields 
rely  heavily  on  the  use  of  the  electromagnetic  spectrum  (EMS),  whether  for 
surveillance  and  target  acquisition,  passage  of  information,  processing  of 
information,  or  destruction  of  enemy  forces  (Frater,  &  Ryan,  2001).  Figure  1 
shows  all  the  necessary  capabilities  for  10  and  their  roles. 
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Figure  1.  Information  Operations  Integration  (From:  JP  3.13) 

B.  ELECTRONIC  WARFARE  AND  ITS  COMPONENTS 

Electronic  Warfare  is  one  of  the  five  core  capabilities  of  10.  EW  dominates 
the  EMS  for  the  benefit  of  friendly  forces.  The  EMS  is  shown  in  Figure  2  and 
Figure  3.  EW  can  be  defined  as  "...  any  military  action  involving  the  use  of 
electromagnetic  (EM)  and  directed  energy  to  control  the  EM  spectrum  or  to 
attack  the  adversary.”  (JP  3.13,  2006)  The  means  to  conduct  EW  targeting  are 
typically  technological  in  nature,  and  the  immediate  targets  are  normally 
technological,  but  the  ultimate  target  is  the  commander’s  decision-making 
capability.  If  EW  is  successful,  the  friendly  commander  can  make  good  decisions 
based  on  the  information  coming  from  technological  devices;  on  the  other  hand, 
the  adversary  commander  cannot  make  good  decisions  because  he/she  cannot 
access  good  and  useable  information  as  much  as  he/she  needs  to. 
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THE  ELECTROMAGNETIC  SPECTRUM 


Figure  2.  Electromagnetic  Spectrum  (From  JP  3.1 3.1 ) 


Figure  3.  Electromagnetic  Spectrum  (From:  National  Aeronautics  and  Space 

Administration,  Science  Mission  Directorate) 
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All  components  of  EW  can  be  applied  to  all  kind  of  operations.  During 
peacetime,  military  forces  try  to  use  EW  to  detect  potential  adversary  EMS  usage 
and  gather  intelligence;  during  wartime,  they  use  EW  to  protect  their  own  EMS 
usage  ability  and  prevent  adversary  usage. 

EW  can  be  divided  into  two  parts:  communications  EW  and  non¬ 
communications  EW.  Communications  EW  is  mostly  concerned  with 
communication  sources  that  transmit  in  frequency  bands  between  HF  (3-30 
MHz),  VHF  (30-300  MHz),  UHF  (300-3000  MHz),  and  SHF  (3000  MHz  to  30 
GHz).  Non-communications  EW  is  mainly  concerned  with  radar  systems,  of 
which  some  operate  in  the  lower  communications  bands,  but  most  are  located  in 
the  higher  microwave  and  millimeter  wave  frequencies.  Additionally,  research  in 
directed  energy  weapons  and  technology  is  increasing  in  importance. 

As  shown  in  Figure  4,  EW  has  three  doctrinal  subdivisions. 


Figure  4.  Subdivisions  of  EW 


1.  Electronic  Attack 

EA  involves  actions  taken  against  personnel,  equipment  or  facility  to 
degrade  their  use  of  the  EMS  and  combat  abilities. 

Subdivisions  of  EA  are  shown  in  Figure  5. 
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Jamming 


Figure  5.  Subdivisions  of  EA  (From  Frater  &  Ryan,  2001) 

2.  Electronic  Protection 

EP  involves  actions  taken  to  protect  personnel,  equipment  and  facilities 
from  any  effect  of  friendly  or  adversary  EW  activities  that  degrade,  neutralize  or 
destroy  friendly  combat  capabilities. 

Subdivisions  of  EP  are  shown  in  Figure  6. 
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Siting 


Figure  6.  Subdivisions  of  EP  (From  Frater  &  Ryan,  2001) 

3.  Electronic  Warfare  Support 

Since  this  thesis  focuses  on  locating  adversary  emissions  in  the  difficult 
terrain  of  eastern  Turkey,  our  discussion  of  ES  will  be  more  extensive  and 
relevant  to  the  thesis  than  the  other  doctrinal  subdivisions  of  EW.  ES  involves 
actions  taken  to  identify,  intercept  and  locate  intentional  or  unintentional  radiated 
electromagnetic  energy.  The  purpose  of  the  ES  is  target  recognition.  The  main 
functions  of  ES  are  to  produce  intelligence,  to  produce  steerage  for  EA  and  to 
cue  surveillance  and  target  acquisition  resources.  (Frater  &  Ryan,  2001) 

Subdivisions  of  ES  are  shown  in  Figure  7. 
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Search 


Figure  7.  Subdivisions  of  ES  (From  Frater  &  Ryan,  2001) 

ES  differs  from,  but  is  similar  to  Signals  Intelligence  (SIGINT).  ES  is  used 
for  immediate  battlefield  information  and  SIGINT  is  used  for  intelligence.  ES 
supports  near  term  operational  applications  and  SIGINT  supports  long  term 
applications.  The  combat  information  gathered  by  ES  can  be  provided  to 
intelligence  resources  in  addition  to  being  used  operationally.  Combat 
information  does  not  normally  require  the  type  of  deep  analysis  that  SIGINT 
typically  does. 

Previously  shown  in  Figure  7  were  the  subdivisions  of  ES,  which  we  will 
now  discuss  in  more  detail. 

Search:  it  is  necessary  to  search  for  and  identify  the  EM  signal  that  the 
adversary  uses  in  the  EMS  before  it  can  be  examined. 

The  search  systems  act  in  space,  time,  and  frequency.  They  have  to  be 
close  enough  to  the  adversary’s  transmitting  system  to  be  able  to  detect  the 
signal.  They  have  to  be  actively  searching  or  listening  at  the  same  time  that  the 
adversary  system  is  transmitting.  Finally,  they  must  be  listening  to  the  same 
frequencies  used  by  the  transmitter.  This  frequency  requirement  is  generally  met 
with  two  compatible  technology  approaches,  narrowband  receivers  and 
broadband  receivers.  Narrowband  receivers  can  receive  a  single  signal  at  a  time 
and  can  scan  a  desired  bandwidth  sequentially  in  frequency;  on  the  other  hand, 
broadband  receivers  can  monitor  multiple  channels  at  the  same  time. 
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Intercept:  Once  the  signals  of  interest  are  identified  through  the  search 
process  they  have  to  be  analyzed  during  intercept  based  on  their  modulation, 
bandwidth,  amplitude,  frequency  and  other  parameters.  This  process  is  also 
called  monitoring. 

Direction  Finding  (DF):  The  location  of  the  transmitter  is  determined  by 
information  gathered  during  the  search  and  intercept  process.  These  locations 
are  likely  to  be  an  approximation  rather  than  an  exact  location.  DF  was 
historically  based  on  triangulation  where  there  had  to  be  at  least  three  receivers 
around  the  emitter. 

DF  systems  historically  employed  special  antennas,  which  defined  the 
bearing  towards  the  emitter.  When  the  lines  of  bearing  from  each  receiver  are 
drawn  on  a  map  manually  or  automatically,  the  interception  of  the  lines  form  a 
triangle.  The  smaller  the  triangle,  the  better  the  accuracy  of  the  system.  The 
emitter  lies  inside  the  triangle. 

Analysis:  Once  the  signals  are  examined  they  are  analyzed  to  define  the 
adversary’s  electronic  warfare  capabilities.  The  main  purpose  is  to  clarify  the 
battlefield  for  the  commander  from  the  EMS  perspective. 

ES  targets  an  adversary’s  EA,  communication  systems  and  electronic 
systems.  A  typical  target  is  an  adversary’s  communication  systems,  where  the 
information  gathered  is  used  for  operationally  actionable  intelligence  and 
targeting  purposes. 

Because  JP  3.1 3.1  is  the  shaping  document  of  EW  for  the  U.S.  and  many 
of  its  allies,  it  is  a  good  idea  to  supplement  Frater  and  Ryan’s  book  with 
discussions  from  the  joint  publication.  According  to  JP  3.13.1,  EA  has  five 
subdivisions:  Electromagnetic  Jamming,  Electromagnetic  Deception,  Directed 
Energy,  Anti-radiation  Missiles  and  Expendables  (e.g.,  chaff,  flares  and  active 
decoys).  ES  has  three  subdivisions:  Threat  Warning,  Collection  Supporting  EW 
and  Direction  Finding.  EP  has  three  subdivisions:  Spectrum  Management,  EM 
Hardening  and  Emission  Control.  Subdivisions  of  EW  do  not  act  alone,  they 
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interact  with  each  other.  The  interaction  among  the  subdivisions  of  EW  and  their 
subdivisions  are  shown  in  Figure  8. 


Figure  8.  Overview  of  EW  (From  JP  3.13.1 , 2007) 

Electronic  warfare  can  also  be  categorized  by  whether  it  is  considered 
active  or  passive.  ES  tends  to  be  passive,  EA  tends  to  be  active  and  EP  tends  to 
be  both  passive  and  active.  Figure  9  shows  this  relationship.  Active  activities 
require  emission  of  detectable  signals  by  the  party  conducting  EW  that  are 
transmitted  (such  as  in  the  example  of  the  jamming  of  a  radar).  Passive  activities 
do  not  emit  signals,  but  rather  depend  upon  detection  of  signals  emitted  by  a 
targeted  emitter.  Active  activities  can  be  normally  be  implemented  during 
peacetime  only  under  strict  limitations;  on  the  other  hand,  passive  activities  can 
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be  implemented  during  peacetime  with  few  if  any  limitations  (Frater  &  Ryan, 
2001). 


Figure  9.  Categorization  of  EW  based  on  Active  and  Passive  (Frater  &  Ryan,  2001) 

In  this  Chapter  we  have  discussed  the  importance  of  Electronic  Warfare. 
In  the  next  Chapter  we  develop  the  fundamentals  of  geolocation  that  will  be 
important  to  our  study. 
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III.  EMITTER  GEOLOCATION 


A.  INTRODUCTION  TO  EMITTER  GEOLOCATION 

Determining  the  location  of  a  target  emitter  is  one  of  the  fundamental 
operations  of  EW,  and  can  serve  many  useful  purposes.  The  position  of  an 
emitter  may  indicate  the  position  of  the  enemy  forces.  In  addition,  precise 
location  of  the  target  emitter  enables  the  use  of  Global  Positioning  Systems 
(GPS)  based  weapons. 

For  civil  purposes,  knowing  the  positions  of  nodes  in  a  wireless  network 
for  commercial  uses  enables  a  variety  of  functionalities,  such  as  emergency 
services,  identification  and  tracking,  location  dependent  computing,  health 
monitoring  and  geographic  routing  (Xu,  Ma,  &  Law,  2006). 

Passive  location  finding,  in  addition  to  active  location  finding  techniques, 
can  locate  stationary  and  moving  targets  by  measuring  the  electromagnetic 
radiation  emitted  by  a  target  with  the  added  benefit  of  not  having  to  radiate 
electromagnetic  energy  to  locate  it.  Passive  and  active  location  finding 
technologies  play  important  roles  in  navigation,  aviation,  aerospace  and 
electronic  warfare  (Yan-Ping,  Feng-Xun,  &  Yan-Qui,  2010). 

The  purpose  of  direction  finding  (DF)  is  to  estimate  or  fix  the  position  of 
selected  emitters.  This  position  is  not  certain  because  all  the  measurements 
include  some  sort  of  error,  and  the  entire  system  includes  the  noise  that  is  found 
in  all  communications  systems. 

Several  techniques  can  be  used  to  calculate  the  position  of  a  target 
emitter.  These  techniques  are  based  on  different  types  of  information  acquired 
from  the  received  signal  to  calculate  a  position  fix  (PF). 

The  azimuth  angle  of  arrival  of  a  signal,  or  so  called  line  of  bearing  (LOB), 

is  the  most  commonly  used  technique  for  calculating  a  PF.  Two  or  more  LOBs 

are  used  to  determine  a  position  in  two  dimensions  (2D).  These  LOBs  are 

assumed  to  be  measured  at  the  same  time  on  the  same  target.  These  LOBs  may 
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intersect  as  illustrated  in  Figure  10.  This  technique  is  called  triangulation  (Poisel, 
2005). 

Since  two  LOBs  intersect  at  a  point,  the  information  to  fix  the  position  of 
the  emitter  is  not  accurate  due  to  measurement  and  propagation  errors.  For 
triangulation,  at  least  three  receivers  are  located  on  a  baseline  as  illustrated  in 
Figure  10.  Each  DF  receiver  has  special  antennas.  These  antennas  are  used  to 
measure  the  bearings.  These  bearings  are  plotted  on  a  map  either  manually  or 
automatically.  Intersection  of  these  bearings  (lines)  forms  a  triangle  and  the 
possible  location  of  the  target  is  calculated  to  be  at  the  middle  of  the  triangle.  The 
size  of  the  triangle  depends  on  the  accuracy  of  measurements.  The  smaller  the 
triangle,  the  better  the  accuracy  (Frater  &  Ryan,  2001). 


Receiver  2 


Figure  10.  Intersection  of  measured  LOBs 

Another  emitter  location  technique  is  to  measure  the  time  of  arrival  (TOA) 

of  a  signal  at  several  sensors.  The  TOAs  can  be  used  to  calculate  the  position 

directly,  but  typically  they  are  sent  to  a  central  node  where  the  time  difference  of 
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arrival  (TDOA)  is  calculated  from  every  pair  of  TOA.  Then  the  range  differences 
between  sensors  and  the  target  are  calculated.  These  range  differences  are 
related  to  the  TDOAs  by  the  speed  of  propagation  in  the  medium.  In  the  air,  this 
is  assumed  to  be  the  speed  of  light. 

The  TDOA  technique  generates  quadratic  lines  of  position  (LOP).  All  the 
LOPs  are  subject  to  measurement  errors  and  noise.  The  intersection  of  these 
LOPs  is  used  to  define  the  emitter  position  (Poisel,  2005). 

The  next  step  of  DF  is  geolocation.  Geolocation  is  closely  related  to  DF 
and  triangulation,  but  it  is  more  realistic  and  distinguished  from  DF  by 
determining  a  meaningful  location  rather  than  just  a  set  of  geographic 
coordinates. 

“Emitter  geolocation  has  two  components.  One  is  measurement,  or  choice 
of  sensors,  and  the  other  is  estimation/information  fusion,  or  processing  of 
measurements  provided  by  the  sensors.”  (Musicki  &  Koch,  2008) 

Geolocation  of  a  source  has  a  wide  variety  of  applications,  such  as 
location  of  radar  sites.  Localization  of  a  interference  source  in  satellite 
communication  systems  is  another  example  of  geolocation  (Sathyan,  Kriburajan, 
&  Sinha,  2004). 

Geolocation  is  based  on  techniques  that  rely  on  frequency,  time,  or  spatial 
information,  or  a  combination  of  these.  Common  methodologies  use  angle  of 
arrival  (AOA),  TOA,  TDOA,  or  differential  doppler  (DD),  also  called  FDOA  (Ho  & 
Chan,  1993)  (Loomis,  2007). 

B.  BEARING  ESTIMATION 

Several  techniques  that  can  be  used  to  determine  the  LOPs  are  discussed 
in  this  section.  The  phase  or  time  difference  can  be  measured  between  the 
signals.  The  amplitude  difference  between  two  signals  can  also  be  measured. 
Frequency  differences  measured  can  be  used  to  determine  the  bearing,  if  one  or 


23 


more  of  the  receivers  are  moving  relative  to  the  other  or  to  the  target  (Poisel, 
2005). 

LOP  systems  do  not  operate  at  the  frequencies  of  the  signals.  The 
frequency  is  typically  converted  to  an  Intermediate  Frequency  (IF)  and  the 
phase/time  measurements  are  made  on  this  converted  signal  (Poisel,  2005). 

1.  Circular  Antenna  Array 

One  of  the  common  types  of  antenna  arrays  for  bearing  determination  is  a 
circular  array.  An  example  of  a  four-element  array  is  illustrated  in  Figure  11. 
Other  forms  of  this  antenna  may  include  more  or  fewer  elements.  A  sense 
antenna  is  included  with  a  circular  array.  The  sense  antenna  is  used  to  remove 
ambiguities  (Joong,  Chul-Gu,  &  Gyu,  2004). 

Let  R  be  the  radius  of  the  circular  array.  R  is  the  length  of  an  “arm”  of  the 
array  measured  from  the  center  to  any  one  of  the  antenna  elements.  The  2R/A  is 
referred  to  as  the  aperture  of  a  circular  array  (Poisel,  2005). 


Figure  1 1 .  Four  Element  Circular  Antenna  Array  (From  Poisel,  2005) 
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An  incoming  signal  s(t)  and  its  accompanied  E  field  with  a  magnitude  E  is 
illustrated  in  Figure  12.  The  vertical  component  of  this  E  field  is  the  only  portion 
that  affects  a  vertically  oriented  antenna.  Thus,  Evert  is  given  by  Ecosa  and  the 
corresponding  signal  amplitude  must  be  adjusted  by  this  factor  (Poisel,  2005). 


2.  Interferometry 

One  of  the  techniques  for  measuring  the  AOA  of  a  coming  signal,  and 
determining  its  LOP,  is  interferometry.  In  interferometry  the  phase  difference  or 
time  difference  between  two  antennas  is  measured  directly.  The  distinction 
between  this  technique  and  other  techniques  lies  in  what  is  measured. 
Interferometry  requires  highly  accurate  measurements.  In  order  to  achieve  this 
performance  accuracy  it  is  necessary  to  space  at  least  two  of  the  antennas  such 
that  the  range  of  possible  phase  difference  can  exceed  2tt  radsian  (Vaccaro, 
1993). 
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Phase  comparison  DF  systems  consist  of  several  antenna  elements  which 
are  arranged  in  a  particular  geometric  configuration.  The  number  of  the  elements 
and  the  arrangement  depends  on  the  DOAs  of  interest  and  the  method  used  to 
process  the  signals.  The  minimum  number  of  elements  is  two.  The  determination 
of  DOA  is  performed  by  direct  phase  comparison  of  the  received  signals  from  the 
different  antenna  elements. 

There  are  two  types  of  interferometers.  The  first  type  measures  the  phase 
difference  between  the  two  antennas  and  calculates  the  AOA  from  that 
measurement.  These  interferometers  are  called  phase  interferometers.  The 
second  type  measures  the  time  of  arrival  differences  to  the  two  antennas  and 
calculates  the  AOA  from  these  measurements.  These  types  of  interferometers 
are  called  active  time  interferometers.  For  this  second  type  of  interferometer 
there  must  be  some  sort  of  time  mark  to  be  able  to  measure  the  time  difference. 
If  there  is  no  time  mark  then  the  signal  reaching  these  two  antennas  must  be 
correlated.  Radar  pulses  provide  a  convenient  time  mark  in  their  leading  edge. 

Interferometer  systems  operate  well  over  limited  frequency  ranges.  They 
provide  the  capability  of  receiving  signals  with  good  accuracy  (Lipsky,  1987). 

C.  QUADRATIC  POSITION  FIXING  TECHNIQUES 

Techniques  for  determining  a  position  fix  based  on  Time  of  Arrival  (TOA), 
Time  Difference  of  Arrival  (TDOA)  and  Frequency  Difference  of  Arrival  (FDOA) 
(also  known  as  Differential  Doppler  or  Differential  Frequency)  techniques  are 
explained  in  this  section. 

When  TOA  or  TDOA  is  measured  at  two  or  more  widely  separated 
receivers,  quadratic  LOPs  are  determined.  The  intersection  of  these  LOP  curves 
is  taken  as  the  estimated  location  of  the  emitter. 

The  receivers  and  the  emitter  might  be  stationary  for  TDOA  and  TOA  but 
for  FDOA  either  the  receiver  or  the  emitter  must  be  moving  in  order  to  produce  a 
frequency  difference  induced  by  movement,  necessary  for  measurement. 
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The  main  advantages  of  using  these  techniques  are: 

•  Most  of  the  time  only  one  antenna  is  used  per  receiver  instead  of 
two  or  more  antennas  as  in  interferometric  techniques. 

•  Higher  precision  and  more  accuracy  can  be  obtained  with  these 
techniques  (Poisel.  2005). 

On  the  other  hand,  there  are  two  main  disadvantages: 

•  Preprocessed  data  samples  are  required  for  non-pulsed 
modulation,  like  Continuous  Wave  (CW)  signals  such  as  using 
Frequency  Modulation  or  Amplitude  Modulation  (Poisel,  2005).  For 
pulsed  conventional  radar  signals  it  is  easy  to  define  the  beginning 
of  the  signal,  but  for  CW  signals  more  complicated  methods  like 
correlation  should  be  used. 

•  Also,  it  is  difficult  to  measure  frequency  of  pulse-type  signals  to  the 
level  of  accuracy  needed  to  do  FDOA  because  the  frequency 
resolution  is  equal  to  1/7",  where  Tis  the  pulse  duration. 

This  following  section  begins  with  a  presentation  of  the  TDOA  technique, 
followed  by  a  discussion  of  FDOA,  and  concludes  with  a  discussion  of  TOA. 

1.  Time  Difference  of  Arrival  (TDOA) 

The  TDOA  involves  the  measurement  of  the  TOA  of  the  received  signal. 
The  TDOA  technique  needs  two  or  more  geographically  separated  sensors 
synchronized  with  each  other  to  be  able  to  find  the  location  of  the  emitter.  If  only 
two  receivers  are  used,  this  will  normally  result  in  an  ambiguous  solution  of  two 
locations,  and  a  third  receiver  is  necessary  to  resolve  this  ambiguity.  These 
receivers  can  be  either  on  the  ground  or  on  an  airborne  platform. 

The  geometry  of  TDOA  is  shown  in  Figure  13.  Only  two  receivers  are 
shown  for  simplicity  and  will  be  used  for  calculations.  This  geometry  is  for  2D  but 
if  the  receivers  or  the  emitter  is  elevated  it  becomes  3D.  For  both  the  scenarios,  r 
is  the  slant  range  where  o  is  the  range  between  receiver  1  and  emitter  and  r2  is 
the  range  between  receiver  2  and  emitter. 
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Figure  1 3.  Geometry  of  TDOA 


The  geometry  involves  an  emitter  at  position  x  =  (xT,yT)  and  two  receivers 

at  positions  (-^,0)  and  (x2,0),  At  time  t  =  0,  a  measurement  of  the  TDOA  is 

made  between  the  arrival  of  the  same  pulse  from  the  emitter  arriving  at  the  two 
receivers. 

The  distances  between  receivers  and  the  emitter  can  be  calculated  as 
follows: 


r,  =  ctt  i  =  1,2 


(3-1) 


where; 

c  is  the  speed  of  light, 

tj  is  the  time  between  when  the  signal  leaves  the  emitter  and 
when  it  arrives  the  receiver. 
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TDOA  is  the  time  difference  between  when  the  signal  arrives  at  one 
receiving  site  and  at  the  other,  and  represented  as  r 

r  =  4-*i=- =  (3.2) 

c  c  c 

r,.J(xT  -  Xf  +  yT2  '=1.2  (3.3) 

re  =  5=  ^(xr - x,)2  +  yr2  - ^(xr - x2)2  +  yr2  (3.4) 

Squaring  both  sides  of  the  previous  equation  yields:  (the  receivers  are  at 
the  same  distance  from  the  origin) 

52  =  2xr2  +  2y/  +  2s2  -  2,|(xr  +  s)2  +  yr2  ]  [( xT  -  s)2  +  yr2  ]  (3.5) 

where 

X,  =  -s  &  x2  =  s 

With  some  more  algebra,  this  expression  becomes 

-  4s2S2  =  (4S2  - 1 6 s2)xt2  +  4 82yT2  (3.6) 

After  further  manipulation,  this  becomes: 
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(3.7) 


xT 


Yt 


82! 4  (4s  -8  )/4 


This  is  familiar  equation  for  a  hyperbola,  which  has  x-axis  intercept  of  (0, 
0)  and  is  asymptotic  to  the  lines  (Loomis,  2007) 


y  =  ± 


^4s2  -  8 2 
8 


x 


(3-8) 


The  curve  defined  by  this  expression  is  a  hyperbola.  Several  of  the 
hyperbolic  curves  are  shown  in  Figure  14. 


Figure  14.  Hyperbolic  TDOA  Curves  (From  Loomis,  2007) 
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It  is  clear  that  two  receivers  cannot  find  the  location  of  the  emitter  since 
these  hyperbolas  covers  a  wide  variety  of  locations.  At  least  three  receivers  are 
necessary  to  find  the  geolocation  of  the  emitter  in  two  dimensions  or  on  the 
earth’s  surface. 

The  following  derivation  of  the  solution  for  the  location  of  the  emitter  as  a 
mathematical  model  of  TDOA  is  quoted  from  Poisel’s  book,  Introduction  to 
Communication  Electronic  Warfare  Systems,  Chapter  12. 

Suppose  there  are  S  receivers  available  to  receive  and  compute  the 
location  of  the  emitter,  then  the  Equation  (3.1)  becomes  for  all  pairs  of  sensors 
(Poisel,  2002). 


di-dj  =  c(ti-tj)  =  ctj 


/,y  =  0,1,...S-1  (3.9) 


According  to  Poisel’s  approach,  let’s  assume  that  all  of  the  arrival  times 
are  compared  with  the  arrival  time  at  a  sensor  located  at  coordinates  (0,0)  as 
shown  in  Figure  15. 
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Figure  1 5.  Sensor  Grid  and  target  in  Two  Dimensions 


The  time  differences  of  arbitrary  (i,j)  are  not  used,  just  (/, 0)  is  used  for  all 
/.The  Equation  (3.9)  in  this  case  reduces  to 


d,  =  IK  rT  1 1 =  c(t,,  0  +  ^ 


(3.10) 


When  putting  the  locations  of  the  emitter  and  the  receivers  in  to  Equation 
(3.10),  it  becomes 


[*/  y,  zi] 


X-r 


yT 


+ 


ctiflyjxT2  +  yT2  +  ZT2 


(3.11) 


=4cv4(x'2+y'2+z'2) 
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Putting  Equation  (3.11)  into  matrix  form 


PiXT + cy  xT  i =“Cf/i0 


2  1ll 
+  -  P; 
2 11 


(3.12) 


where 


p,  is  the  position  vector  of  receiver  / 
Expanding  this  result  for  all  /  =  1  ...S  receivers  yields 


PxT  +c||xT||  =  d 


(3.13) 


where 


A  *i  /.  zi  A 


*s  i  ys  i  z 


s-1  J 


XT  = 


xT 


1,0 


t  = 


s-1,0 


c  =  ct 
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d  =  ^ diag{PPT  -  c2ttr) 


Let  a  =  ||xT||  and  Q  =  (PPr)  1  ,  then  the  expression  becomes 


(crQc  - 1 )  a2  -  2drQca  +  drQd  =  0 


(3.14) 


which  is  a  quadratic  equation  that  can  easily  be  solved  for  a  which  is  the 
range  of  the  target  from  the  origin.  Substituting  this  range  back  into  Equation 
(3.13)  will  solve  the  problem  for  the  target  location. 


PxT  +  ac  =  d 

xT  =P  ^(d-ac) 


(3.15) 


If  S>4,  this  becomes  an  over-determined  system  of  equations.  Instead  of 
the  inverse,  then,  the  pseudo  inverse  can  be  used  and  is  shown  in  Equation 
(3.16) 


pt  =(P7P)-1pr 


(3.16) 


The  pseudo  inverse  solves  for  the  emitter  position  in  the  minimum  least 
squares  error  sense  (Poisel,  2002).  It  is  very  important  for  this  thesis  and  for 
many  of  TDOA  solutions.  The  pseudo  inverse  is  used  in  Chapter  4  for  simulation 
which  is  a  developed  version  of  Ezzat’s  Closed-Form  Geolocation  Solution. 
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2.  Frequency  Difference  of  Arrival-Differential  Doppler  (FDOA-DD) 

The  signal  emitted  by  a  target  of  interest  which  is  moving  produces  an 
effect  called  Doppler  shift.  The  Doppler  shift  is  related  to  the  direction  of  the 
movement  relative  to  the  receiver  and  shows  itself  as  a  frequency  difference.  The 
movement  can  be  at  the  target  or  at  the  receiver.  Each  of  these  movements 
produces  the  same  frequency  difference  effect.  When  the  frequency  difference 
between  the  target  and  the  receivers  of  two  (the  more  the  receivers  the  better  the 
calculation)  is  measured,  these  measurements  can  be  used  to  calculate  the 
geolocation  of  the  target.  This  type  of  geolocation  technique  is  called  frequency 
difference  of  arrival  (FDOA)  or  differential  Doppler  (DD). 

Since  this  thesis  is  interested  in  stationary  receivers  and  FDOA  technique 
requires  a  moving  receiver  or  a  moving  target,  FDOA  is  not  examined  in  detail. 

D.  CLOSED-FORM  SOLUTION  OF  HYPERBOLIC  GEOLOCATION 

EQUATIONS 

3D  geolocation  is  important  for  the  Turkish  Army  because  the  terrain 
where  most  of  the  anti-terrorist  operations  take  place  is  very  mountainous, 
especially  the  southeast  part  of  Turkey.  It  is  not  possible  to  locate  or  position  the 
receivers  at  the  same  elevation  with  the  transmitter.  It  is  inevitable  to  be 
presented  with  a  3D  problem.  Ezzat’s  approach  gives  a  unique  solution  to  a  3D 
problem  from  the  TDOA  perspective.  Because  of  that  reason,  Ezzat’s  approach 
is  used  for  developing  and  analyzing  the  scenarios  in  Chapter  4. 

The  following  approach  is  paraphrased  from  Ezzat’s  article. 

The  approach  that  is  explained  here  does  not  depend  upon  range  data. 
Range  data  is  derived  by  multiplying  the  time  that  the  signal  propagated  in  air 
and  the  speed  of  the  signal,  the  speed  of  light.  It  is  necessary  to  know  the  time  of 
transmission  and  the  time  of  arrival  to  get  the  amount  of  time  travelled  from  the 
transmitter  to  the  receiver.  For  the  most  part,  it  is  impossible  to  know  the  time  of 
transmission,  which  can  be  defined  as  f0.  Ezzat  indicate  in  his  paper  that  this  is 

the  only  technique  that  can  be  used  without  the  range  data.  He  mentions  that 
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without  t0,  it  is  possible  to  find  the  location  of  the  emitter.  He  also  indicates  that 
the  previous  solutions  are  good  for  a  noise  free  environment. 

The  closed-form  solution  presented  does  not  require  the  calculation  of 
range  data  and  does  not  depend  on  the  availability  of  any  information  other  than 
the  times  of  arrival. 

The  basic  form  of  time  of  arrival  equation  is  as  follows 


''=^t 


(3.17) 


where 

f,  is  the  time  of  arrival  at  receiver  /, 

D,  is  the  distance  between  the  emitter  and  the  receiver, 
c  is  the  speed  of  light. 
t0  is  the  time  of  the  transmission 

Two  or  more  receivers  are  needed  to  be  able  to  calculate  the  TDOA  as 
mentioned  before.  When  this  condition  is  satisfied,  it  is  possible  to  eliminate  t0 
from  any  pair  of  two  equations,  which  results  in 


=  (3.18) 

This  is  the  TDOA  equation.  This  equation  yields  a  3D  hyperboloid  as 
mentioned  before.  When  the  emitter  coordinates,  x0,  y0,  and  z0,  are  plugged 

into  this  equation, (3. 18)  can  be  written  as 
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>/(* 2  ~xo)2+ (y2  -y0)2+(z 2  -  zo  )2  -  >/(*i  -  *0  )2  + (/i  -  y0  )2 + (zi  -  zo  )2 


=  c(f2-0 


(3.19) 


where 


^ ,  y1;  z,  and  x2,  y2,  z2  are  the  coordinates  of  the  receiving 
antennas  1  and  2. 

Having  three  more  receiving  antennas  yields  three  additional  times  of 
arrival,  which  produce  an  additional  three  equations  like  Equation  (3.19),  it  is 
possible  to  solve  for  the  emitter  coordinates  x0,  y0 ,  z0 

There  is  an  effect  of  path  delay  on  (3.18)  and  hence(3.19).  Ezzat  says  that 
a  propagation  mode  between  any  two  points  in  which  the  path  of  the  signal  is  not 
a  straight  line  will  be  mathematically  equivalent  to  propagation  along  a  straight 
line  but  with  a  velocity  that  is  less  than  c,  as  the  time  of  arrival  is  important.  The 
TDOA  equation  for  the  case  in  which  the  path  of  the  signal  is  nonlinear,  it  will  be 
written  as  follows  after  adding  the  path  delay: 


D, 


a2c  a^c 


1(4  _9i) 

C  a2 


(3.20) 


where 

a1  and  a2  are  path  delay  coefficients  which  are  less  than  or  equal 

to  1 . 

a  is  one  when  there  is  no  path  effect  on  the  propagating  signal  like  in  the 
case  where  the  signal  propagates  through  air.  For  this  thesis,  all  the  path  delay 
coefficients  are  assumed  to  be  one. 
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The  solution  that  is  used  for  this  thesis  is  derived  from  Ezzat’s  article. 
First,  I  will  explain  his  approach.  Then,  because  he  did  not  discuss  the  effects  of 
noise,  I  will  add  noise  to  his  approach. 

These  are  the  steps  to  transform  the  hyperbolic  equations  into  a  set  of 
vector  equations.  In  Equation  (3.18),  we  first  note  that  a  distance  dj  can  be 
written  as  the  norm  of  a  vector. 


cHPi-Po 


(3.21) 


where 

Pi  =(xl,yi,zi) 

and 


is  the  position  vector  of  the  receiving  antenna 


Po  =  (Wo’zo)  is  the  position  vector  of  the  emitter. 


Equation  (3.20)  represents  the  difference  between  the  two.  We  modify 
Equation  (3.20)  and  write  it  as  a  difference  of  squares 


1  rJ  ^  ri  ^ 
(t2-t0)2-(t,-t0)2  =  M\-%) 


an 


a , 


(3.22) 


From  Equation  (3.21)  and  Equation  (3.22)  we  have 
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(3.23) 


=  c2 (42  +  fo2  - 2M,  -  f.2  -  'o2  +  2 Vo) 

a2  v  7 

=  c2(t22-t?)-2t0c2{t2-ti) 

Pi  -  p0 12  can  be  written  as 


Pi-Pof  =  |Pi|2  +  |Po|2  -  2PiTPo 


(3.24) 


where 


represents  the  transpose  of  vector  Pj. 

If  we  put  Equation  (3.24)  into  Equation  (3.23)  it  gives 


p2  r + ip0  r  ■  2p2tpo  iPi  r + ip0  r  -  2pitp0 


a0 


a . 


c2(t22-t2)-2t0c2(t2-t,)  (3.25) 


The  two  coefficients  a1  and  a2  are  both  very  close  to  one,  and  as  a  result 

the  difference  |p2 12  /a22-  |p1 12  la 2  is  negligible  by  comparison  with  the  other 
quantities.  Then  Equation  (3.25)  becomes 


P2  Pi 


an 


a , 


2Po  ~  ^V)  =  c"  ( 42  -  0 )  ~  2f0< c2 ( 4  -  0 

a2  a,  v  ’ 


(3.26) 
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In  Equation  (3.26)  there  are  two  unknowns,  t0  and  p0 .  To  solve  the 
coordinates  of  the  receiver  (p0)  we  need  three  linearly  independent  equations 

like  Equation  (3.51).  Equation  (3.26)  shows  the  relationship  between  receivers  1 
and  2.  The  other  two  linearly  independent  equations  can  be  between  1  and  3  and 
1  and  4.  This  pairing  of  receivers  makes  the  equations  linearly  independent.  As 
can  be  seen  for  three  equations  at  least  four  receivers  are  required.  For  the 
closed-form  approach,  it  is  necessary  to  get  rid  of  t0.  To  do  this  these  linearly 

independent  equations  are  going  to  be  coupled  to  subtract  from  each  other.  For 
that  purpose  this  solution  needs  five  receivers  to  work. 

Equation  (3.26)  is  divided  by  {t2  -t,)  ,  then  it  becomes 


1 

( \  I2 

lp2l 

1  I2  ^ 

|Pi| 

o 

a 

CM 
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(n  T 

P2 

_P^ 

(' 2-4) 

a22 
v  ^ 

a.2 

1  J 

1  (4-01 

V  a2 

«i2v 

=  c2{t2  +  t,)-2t0c2 


(3.27) 


By  similar  steps  Equation  (3.27)  can  be  expressed  for  receiver  pairs  3-1 
and  4-1 


1 

(  1  I2 

|P3| 

1  I2  ^ 

|Pi| 

o 

a 

CM 

1 

(  _  T 

P3 

_pT] 

(h-t,) 

a32 
\  J 
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1  J 

i  te-o1 

V  a3 

a2; 

Equation  (3.28)  is  for  3-1  receiver  pair. 

1 

(\  I2  1  |2  X 

|P4|  |Pl| 

2p0 

f  _  T  nT\ 

P4  Pi 

a4  a ,2 

v  4  1  J 

('4-0 

2  2 
ai  ) 

c2{t3  +  t,)-2t0c2 


c2{t4  +  t,)-2t0c2 


(3.28) 


(3.29) 


40 


Equation  (3.29)  is  for  4-1  receiver  pair. 

The  next  step  is  to  eliminate  t0  from  Equations  (3.27),  (3.28)  and  (3.29). 
(3.27)  and  (3.28)  are  merged  in  the  components  of  p0to  get  rid  of  t0.  This  results 
in 


1 

f  |  |2 

|P2| 

1  I2  Y 

|Pl| 

1 

(  \  |2 
|P3| 

i  i2  A 

|Pi| 

2 

^  a2 

a,2 

1  J 
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^  a3 

a,2 

1  J 

(3.30) 
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Similarly  Equations  (3.27)  and  (3.29)  yield 
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(3.31) 


The  last  two  equations  are  linearly  independent  in  the  components  of  p0 . 

As  mentioned  before  at  least  three  independent  equations  are  required.  For  this 
purpose  there  must  be  another  receiver:  receiver  5.  If  we  follow  the  same  steps 
that  we  followed  to  reach  Equations  (3.30)  and  (3.31),  we  can  have  the  following 
equation  for  receiver-emitter  pair  5-1 . 
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«12  J. 

(3.32) 


Equations  (3.30),  (3.31)  and  (3.32)  can  be  written  in  alternative  algebraic 

form 


3,1*0  ^12/0  ^13^0 

^21*0  ^22X0  ^23  Z0  —  ^2 


^31^0  ^32  34)  ^33^0  ^3 


(3.33) 


where 
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42 
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6,= 


2  ,  .,2  ,  „2  A 


x2  +y2  +z2  ^  +y,  +z, 


a0 


a. 


(  „  2 


( *3  0 


2  ,  .,2  ,  ,2\ 


*3  +  y3  +  zs  *r  +  rr  +  zi 


a,, 


a . 


(3.35) 


The  other  constants  can  be  derived  from  the  same  Equations  (3.31)  and 

(3.32). 

Equation  (3.33)  can  be  expressed  as 


AX  =  B 


(3.36) 


and  can  be  solved  as  in  least  square  estimation  way  as  shown  in  Section 
E. 


X  =  ( AT  *  A)'1  AT  *  B 


(3.37) 


E.  LEAST-SQUARE  ESTIMATION 

For  the  least-square  estimation,  we  use  Barkat’s  book,  Signal  Detection 
and  Estimation,  2005. 

In  the  least-square  estimation,  the  criterion  is  to  minimize  the  squared 
difference  between  the  given  data  (signal  plus  noise)  and  the  assumed  signal 
data. 

This  development  applies  to  a  linearized  version  of  non-linear  equations 
relating  an  /^-dimensional  measurement  vector  and  the  /W-dimensional  vector  to 
be  estimated.  This  linearized  matrix  equation  represents  the  difference  between 
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the  actual  measurement  vector  and  the  measurement  that  would  be  obtained  if 
the  vector  to  be  estimated  has  the  estimated  value. 

Suppose  we  want  to  estimate  M  parameters,  denoting  the  /W-dimensional 
vector  0 ,  from  the  K  measurements,  denoting  the  K'-dimensional  vector  y  with  K 
>  M.  The  relation  between  the  parameters  0  and  the  observed  data  y  is  given 
by  the  linear  model 

Y  =  H0  +  N 

where 

h  is  a  known  (K  x  M)  matrix 

n  is  the  unknown  (K  x  1)  error  vector  that  occurs  in  the 
measurement  of  0 . 

The  least-square  estimator  of  0  chooses  the  values  that  make  X  =  H0 
closest  to  the  observed  data  Y-  Hence,  we  minimize 


J(0)  =  £( Yk  -  Xk  f  =  ( Y  -  H0)r  (Y  - H0)  =  YYT  -  2YTH0  +  0THTH0  (3.38) 


k= 1 


Note  that  YT  -H0  is  a  scalar.  Taking  the  first-order  partial  derivative  of  the 
cost  function  J(9)  with  respect  to  0  and  setting  it  equal  to  zero,  we  obtain  the  set 
of  linear  equations 


=  -2HtY  +  2HTH0  =  0 
dQ 


(3.39) 


and  LSE  can  be  found  to  be 
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0|S  =(HTH)'1HTY 


(3.40) 


We  observe  that  the  error  in  the  estimator  0|S  is  a  linear  function  of  the 
measured  errors  n,  since 


Me-0,s  =  e- (hth y hty  =  0 - (hth)_1  ht [ho + n]  =  -<hth)_1htn  (3.41 ) 

1.  Recursive  Least-Square  Estimator 

In  real  time  estimation  problems,  it  is  necessary  to  write  the  estimator  0  in 

a  recursive  form  for  better  efficiency.  Consider  a  situation  where  an  estimate  0  is 
determined  based  on  some  data  YK.  If  new  data  YK+1  are  to  be  processed  after 

having  determined  an  estimate  based  on  the  data  YK,  it  is  best  to  use  the  old 
solution  along  with  the  new  data  to  determine  the  new  least-square  estimator 
(Barkat,  2005).  “It  is  clear  that  discarding  the  estimate  based  on  the  data  YK  and 

restarting  the  computation  for  a  solution  is  inefficient.  This  procedure  of 
determining  the  least-square  estimate  from  an  estimate  based  on  YK  and  the 
new  data  YK+1  is  referred  to  as  sequential  least-square  estimation,  or  more 
commonly  recursive  least-square  (RLS)  estimation.”  (Barkat,  2005) 

Consider  the  problem  of  estimating  9  from  the  data  vectors  zM  given  by 
the  linear  model 


zM=ly?  +  uM 


(3.42) 


where 


45 


(3.43) 


z  =[Y  Y  Y  V 

is  an  (MK+1)  collections  of  vectors  Y,,Y2,...,YM  since  each  vector 
Yk,  k  =  1,  2,  Mis  a  (K+ 1)  vector 

u  „  =  [N,N2...NuY  (3.44) 

is  an  (MK+1)  error  vector,  and 

hM  =[^i  h2...hMf  (3.45) 

is  an  (MK*n)  mapping  matrix  relating  ZM  to  the  (n*1)  parameter  vector  0to 
be  estimated. 

It  can  be  shown  that  the  RLS  estimator  is  given  by 

-  ®m-i  +Vm[um  (3.46) 

where 

VM  =  CyuhMTRMM  1  (3.47) 

c  is  the  error  covariance  matrix  given  by 
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(3.48) 
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In  three  dimensions,  when  the  covariances  are  fixed  this  matrix  is 


Cuu  - 


Pxyaxay 

Pxz°x°z 


P xy® x^y 
PyzVyVz 


Pxz°xaz 

PyzVyVz 


where 

ax  is  the  variance  of  x 
<ry  is  the  variance  of  y 
<jz  is  the  variance  of  z 

pxy  is  the  correlation  coefficient  between  x  and  y 

pxz  is  the  correlation  coefficient  between  x  and  z 

pyz  is  the  correlation  coefficient  between  y  and  z 

If  it  is  assumed  that  x,  y,  and  z  are  uncorrelated  and 
pxy,pxz,pyz  =  0 .  In  this  case  Equation  (3.49)  becomes  (Poisel,  2005) 


C 


uu 


erxz  0  0 

0  <  0 

0  0  az2 


(3.49) 


so  that 


(3.50) 
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F. 


SOURCES  OF  ERRORS  AND  MEASUREMENT  ACCURACY 


1.  Sources  of  Errors 

Direction  Finding  is  subject  to  number  of  errors.  These  error  sources  are 
listed  and  described  below. 

a.  Equipment  Error 

Modern  DF  equipment  gives  bearing  with  accuracy  of  ±2°.  Hand¬ 
held  tactical  units  have  less  accuracy  of  ±10°  (Frater  &  Ryan,  2001). 

b.  Short-baseline  Error 

If  the  angle  between  bearing  lines  is  less  than  45°,  the  triangle  of 
error  for  triangulation  or  confidence  ellipse  for  hyperbolic  geolocation  systems 
like  TDOA  becomes  significantly  larger  (Frater  &  Ryan,  2001). 

c.  Co-channel  Interference 

Most  tactical  DF  systems  cannot  identify  the  difference  between 
multiple  received  signals.  When  there  is  significant  co-channel  interference  these 
DF  systems  tend  to  give  erroneous  bearing  information  (Frater  &  Ryan,  2001). 

d.  Adjacent  channel  Interference 

“Strong  signals  in  a  channel  adjacent  to  the  one  being  DF’ed  can 
lead  to  an  erroneous  bearing”  (Frater  &  Ryan,  2001). 

e.  Multipath  Error 

In  multipath  error,  two  or  more  signals  arrive  at  the  receiver.  These 
received  signals  originate  from  the  same  source  but  travel  in  a  different  path  due 
to  natural  or  man-made  obstacles  and  reach  the  receiver  at  different  times  due  to 
the  difference  in  the  distance  traveled  over  the  separate  paths. 
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f.  Night  Effect 

Night  effect  is  a  special  case  of  multipath  effect  that  occurs  when 
sky-wave  propagation  occurs  at  night  but  not  during  the  day.  This  type  of  error 
occurs  at  long  distance,  typically  over  HF  communication  channels. 

g.  Coastal  Refraction 

Surface  wave  propagation  that  crosses  a  coastline  at  an  angle 
other  than  a  right  angle  is  subject  to  bending  caused  by  refraction.  This  may  lead 
to  wrong  bearing  determination.  Coastal  refraction  is  usually  significant  at 
frequencies  below  10MHz  (Frater  &  Ryan,  2001). 

h.  Thunderstorms 

Thunderstorms  can  lead  to  a  wrong  bearing  that  points  towards  the 
thunderstorm  rather  than  the  originating  transmitter. 

i.  Rain 

Heavy  rain  may  reduce  received  signal  levels  in  SHF  and  higher 
bands.  This  reduces  the  range  of  DF  systems  (Frater  &  Ryan,  2001). 

2.  Cross-Correlation  TDOA  Estimation  Technique 

The  TDOA  position  fixing  technique  includes  two  phases.  The  first  phase 
is  the  estimation  of  the  TDOAs  of  the  signal  from  a  source,  between  pairs  of 
receivers  through  the  use  of  time  delay  estimation  techniques.  In  the  second 
phase,  the  estimated  TDOAs  are  transformed  into  range  difference 
measurements  between  stations,  resulting  in  a  set  of  hyperbolic  equations 
(Aatique,  1997). 

The  previous  sections  related  to  TDOA  calculations  are  for  the  second 

phase. 

There  are  two  general  methods  for  estimating  the  TDOAs.  The  first  one  is 
to  subtract  TOA  measurements  from  two  stations  to  produce  a  relative  TDOA. 
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The  second  one  is  to  employ  a  cross-correlation  technique,  in  which  the  received 
signal  in  one  station  is  correlated  with  the  received  signal  at  another  station. 

Because  it  is  very  difficult  to  know  the  timing  reference  on  the  source  to  be 
located,  and  because  the  signals  of  interest  for  this  thesis  are  CW,  the  cross 
correlation  technique  is  commonly  used  to  estimate  the  TDOAs.  The  basic  timing 
requirement  for  this  technique  is  to  synchronize  the  receivers.  This  requirement 
is  relatively  easy  compared  to  the  need  to  know  the  originating  transmission  time 
of  the  signal.  Therefore,  we  will  focus  on  the  cross  correlation  TDOA  estimation 
technique. 

Signal,  s(t) ,  emanating  from  a  remote  source  through  a  channel  with 
noise,  the  general  model  for  the  time-delay  estimation  between  received  signals 
at  two  base  stations,  x^and  x2(t) ,  is  given  by  (Knapp  &  Carter,  1976) 


x1(0  =  s(0  +  a?1(0 
x2(t)  =  As(t-r)  +  n2(t) 


(3.51) 


where 

n,(t)  and  n2(t)  are  noises 

t  is  the  TDOA  between  the  receivers 

A  is  the  amplitude  ratio  for  scaling  the  signal 

This  model  assumes  that  s(t),  n^(t)  and  n2(t)  are  real  and  jointly 
stationary  random  processes.  The  signal  s(t )  is  assumed  to  be  un  correlated 
with  noise  n,(t)  and  n2(t) . 

The  cross  correlation  of  this  two  received  signal  is  given  by  (Knapp  et  al, 

1972) 
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^(r)  =  E[x,(Ox2(f-r)] 


(3.52) 


where  e  represents  the  expectation.  Equation  3.51  can  also  be  expressed  as 
(Aatique,  1997) 


=  FP^t)  =  \^x,{t)x2{t-T)dt 


(3.53) 


Because  the  observation  time  cannot  be  infinite  and  can  be  estimated 
from  a  finite  observation  time,  an  estimate  of  the  cross-correlation  is  given  by 


^(r)  =  yjorx1(0x2(f-r)df 


(3.54) 


where  t  represents  the  observation  interval.  The  time  delay  causing  the 
cross  correlation  peak  is  an  estimate  of  TDOA,  r  . 

The  cross  correlation  technique  is  affected  by  many  errors,  which  can  be 
considered  as  a  Gaussian  distribution.  Standard  deviation  for  that  Gaussian 
distribution  is  explain  in  Section  3. 

3.  Standard  Deviation 

When  errors  producing  ambiguity  are  taken  into  consideration,  TDOA 
isochrones  are  no  longer  clear  functions;  they  form  regions  or  areas  where  the 
target  should  be,  as  illustrated  in  Figure  16. 

All  of  the  TDOA  methods  are  subject  to  errors  in  measurement.  The  noise 
and  the  measurement  error  are  the  two  primary  error  sources  (Poisel,  2005).  For 
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the  purpose  of  this  thesis,  errors  discussed  from  this  point  forward  will  be  limited 
to  noise-induced  errors,  and  they  are  discussed  in  terms  of  standard  deviation. 


Figure  16.  TDOA  Isochrones  with  Errors 

The  represented  TDOA  hyperbolic  isochrones  can  be  shown  in  detail  as  in 
Figure  17.  The  standard  deviation  increases  close  to  the  edges  of  the  hyperbolic 
isochrones. 


Figure  17.  TDOA  Isochrones  with  Standard  Deviation 
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The  standard  deviation  discussed  in  this  chapter  result  from  cross¬ 
correlation  measurements. 

The  Cramer-Rao  bound  on  parameter  estimation  is  a  frequently  used 
measure  on  how  well  such  a  parameter  can  be  measured.  The  Caramer-Rao 
bound  for  estimating  the  time  of  arrival  of  a  signal  at  a  receiver  is  given  by  (Stein, 
1981). 


1  1 

IWr 


(3.55) 


where 


B  is  noise  bandwidth  of  receivers 

T  is  integration  time,  which  must  be  less  than  or  equal  to  signal 

duration 

y  is  the  effective  input  SNR  at  two  sensor  sites 

RMS  radian  frequency  is  given  by  p  which  is  the  measure  of  the 
bandwidth  of  the  signal  an  given  by 


p  =  2x 


OT'iTtf  2 


(3.56) 


where 


S(f)  is  the  spectrum  of  the  signal 


Variable  y  is  a  composition  SNR  at  the  two  sensors  and  is  given  by 
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(3.57) 
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where 


yi  and  y.  are  the  Signal  to  Noise  Ratio  (SNR)  at  two  receivers. 
For  low  SNR,  standard  deviation  is  given  by 


1111 


cr 


o7i2  sJTWf0 


i  + 


w2 

12  f: 


(3.58) 


where 

T  is  the  integration  time 

8  is  the  SNR 

W  is  the  bandwidth  and  is  given  by  W  -  f2- 

f0  is  the  center  frequency 

For  high  SNR,  standard  deviation  is  given  by 


cr. 


1  1  1 

4x2T  8  Jf*  _  ft 


(3.59) 


The  standard  deviation  function  for  low  SNR  is  illustrated  in  Figure  18  and 
Figure  19  for  short  integration  time  and  high  integration  time,  and  in  Figure  20 
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and  Figure  21  for  high  SNR  for  the  25MHz  bandwidth  (Poisel,  2005).  The  dotted 
lines  represent  the  function  in  Equation  (3.55)  and  the  solid  lines  represent  the 
function  in  Equation  (3.59). 


Figure  18.  TDOA  Standard  Deviation  for  Low  SNR,  W=25  MHz  and  Long  Initegration 

Time  (From  Poisel,  2005) 
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Figure  19. 


Figure  20. 


TDOA  Standard  Deviation  for  Low  SNR,  W=25  MHz  and  Short  Integration 

Time  (From  Poisel,  2005) 


SNR  (dB) 

TDOA  Standard  Deviation  for  High  SNR,  W=25  MHz  and  Long  Integration 

Time  (From  Poisel,  2005) 
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Figure  21 .  TDOA  Standard  Deviation  for  High  SNR,  W=25  MHz  and  Short  Integration 

Time  (From  Poisel,  2005) 

4.  TDOA  Dilution  of  Precision 

TDOA  measurement  suffers  from  another  type  of  error,  caused  by  the 
long  ranges  from  the  sensor  baseline.  Consider  the  Figure  22  where  the  target  is 
very  far  from  the  baseline  between  the  sensors.  The  hyperbolic  LOPs  are  nearly 
parallel  to  each  other  and  make  the  measurement  more  difficult,  thereby  making 
the  system  vulnerable  to  any  measurement  or  noise  error.  This  is  called 
geolocation  dilution  of  precision  (GDOP)(  Poisel,  2005). 
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Figure  22.  Geometry  of  GDOP  Effect  on  TDOA  (Poise!,  2005) 


The  GDOP  effect  is  illustrated  in  Figure  23  where  the  effect  of  the  distance 
to  the  error  can  be  seen  easily.  N  represents  the  number  of  receivers  in  the 
system. 
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Figure  23.  TDOA  GDOP  from  4  to  10  elements  (From  Bard  &  Ham,  1999) 

5.  Effects  of  Movement  on  TDOA 

When  there  is  a  motion  between  the  receivers  or  the  target  the  error  of 
TDOA  measurements  increases.  This  is  illustrated  by  Chan  and  Ho,  who  called 
this  effect  scale  difference  of  arrival  (SDOA).  Figure  24  shows  the  relationship 
between  the  velocity  and  the  error  caused  by  this  velocity.  As  can  be  seen,  when 
the  velocity  increases  the  mean  square  error  also  increases. 
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Figure  24.  SDOA  Caused  by  the  Relative  Motion  (From  Chan  &  Flo, 2003) 

In  this  Chapter,  we  discussed  the  fundamentals  of  emitter  geolocation, 
bearing  estimation,  quadratic  position  finding  methods,  closed-form  geolocation 
of  emitter  based  on  TDOA,  least-square  estimation  method  and  sources  of  error 
for  geolocation.  In  the  next  chapter,  we  developed  a  Matlab1  simulation  based  on 
Ezzat’s  closed-form  geolocation  technique  and  made  four  type  of  analysis;  the 
effect  of  the  distance  and  noise  on  accuracy,  comparison  between  stationary 
receivers  and  moving  receivers  based  on  their  effects  on  accuracy,  the  optimum 
distribution  of  stationary  receivers  for  best  accuracy. 


1  Matlab  is  a  registered  trademark  of  The  Mathworks,  Inc. 
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IV.  SIMULATION  OF  THE  CLOSED-FORM  GEOLOCATION 

TECHNIQUE 


The  simulation  described  in  this  chapter  focuses  on  the  use  of  time 
difference  of  arrival  (TDOA)  employing  Ezzat’s  proposed  closed-form  solution  for 
the  emitter  position  in  three  dimensions  (Ezzat,  2007).  We  will  analyze  its 
capabilities  in  terms  of  the  effect  of  the  number  of  experiments  on  accuracy,  the 
effect  of  a  noisy  environment  on  accuracy,  the  effect  of  distance  between  emitter 
and  receivers  on  accuracy,  and  the  optimum  positions  of  the  receivers  for  best 
accuracy.  Finally,  we  will  make  a  comparison  between  moving  receivers  and 
stationary  receivers. 

The  assumptions  of  and  the  restrictions  to  the  simulation  are  given  in 
Section  A.  The  development  of  the  simulation  and  the  followed  procedures  are 
explained  in  Section  B.  Simulation  results  are  given  in  Section  C. 

Ezzat’s  closed-form  geolocation  technique  has  been  used  to  develop  this 
simulation.  This  technique  is  explained  in  Section. D  in  detail.  Matlab  has  been 
used  for  simulation,  and  Matlab  code  can  be  found  in  Appendix  A. 

In  short,  simulation  is  done  to  accomplish  these  analyses: 

•  The  effect  of  number  of  experiments  on  the  accuracy  of  the 
estimated  emitter  position, 

•  The  effect  of  the  distance  between  emitter  and  receivers  on  the 
accuracy, 

•  The  effect  of  the  noise  on  the  accuracy, 

•  Comparison  between  moving  receivers  and  the  stationary 
receivers, 

•  The  optimum  geographical  distribution  of  the  receivers  for  best 
accuracy. 

Different  codes  are  developed  for  each  analysis  and  can  be  found  in  the 
Appendices. 
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A. 


ASSUMPTIONS  AND  RESTRICTIONS 


The  following  assumptions  and  restrictions  apply  to  the  development  of 
the  technique. 

1.  Assumptions 

1)  The  sensor  locations  are  known  exactly. 

2)  The  transmitter  location  is  fixed  during  the  period  of  DF  fixing. 

3)  The  sensors  are  properly  located  and  operated. 

4)  The  emitter  is  within  the  range  of  the  receivers. 

5)  TDOA  is  measured  using  some  sort  of  technique  and  known. 

6)  The  error  in  the  noise  is  distributed  as  a  Normal  distribution  with 
zero  mean  and  known  standard  deviation. 

7)  The  noise  for  each  pair  is  independent. 

8)  There  is  no  multipath  effect. 

9)  The  emitter  location  estimate  error  is  distributed  as  a  normal 
distribution. 

Assumption  #1  is  reasonable  based  on  the  fact  that  any  such  position 
errors  can  be  added  to  the  emitter  estimate  uncertainty,  if  they  are  significant. 

Assumption  #2  is  necessary  to  the  analyses  of  the  systems  considered  in 
this  thesis,  and  it  is  reasonable  over  the  period  required  to  obtain  a  single  fix. 
Further,  any  change  in  emitter  position  is  not  significant  in  comparison  to  the 
speed  of  signal  transmission  (speed  of  light)  and  time  to  calculate  solutions. 

Assumption  #3  is  reasonable  in  the  absence  of  contradictory  information. 

Assumption  #4  is  reasonable  based  on  the  fact  that  if  the  emitter  is  out  of 
the  range,  then  the  system  will  not  function. 

Assumption  #5  is  reasonable  because  the  measurement  of  TDOA  is  out  of 
the  interest  area  of  this  thesis. 
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Assumption  #6  is  reasonable  based  on  the  nature  of  the  noise  in  every 
receiver  system. 

Assumption  #7  is  reasonable  based  on  the  fact  that  every  receiver  pair 
has  different  SNR  level,  integration  time  and  bandwidth. 

Assumption  #8  is  reasonable  based  on  the  fact  that  multipath  effect  is  out 
of  the  interest  area  of  this  thesis. 

Assumption  #9  is  usual  when  considering  measurements  which  are 
subject  to  random  measurement  error.  There  are  biases  in  the  measurements 
from  navigation  errors,  errors  in  the  calibration  tables,  interference,  etc.;  these 
biases  may  be  removed.  In  the  absence  of  specific  knowledge  about  these  errors 
the  normal  assumption  is  reasonable.  Biases  will  not  be  treated  in  this  thesis. 

2.  Restrictions 

In  addition  to  the  assumptions  discussed  in  this  section,  this  thesis  does 
not  consider  the  following  effects; 

1)  Geographic  transformation,  map  projection  effects,  and  grid 
reference  system  conversions. 

2)  Propagation  effects. 

3)  Susceptibility  to  deception. 

4)  Special  problems  associated  with  low-probability-of-intercept 
emitters  (low  SNR,  spread-spectrum,  time-frequency  diversity,  frequency  agility, 
etc.). 

5)  Numerical  computation  and  normal  truncation  effects. 

6)  Equipment  errors. 

7)  Interference. 

8)  Night  effect. 

9)  Coastal  refraction. 
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10)  Thunderstorms  and  precipitation. 


B.  DEVELOPMENT  OF  THE  SIMULATION 

Ezzat’s  technique  and  Matlab  code  solve  for  the  position  of  the  emitter 
directly.  The  Matlab  code  includes  multipath  effect  for  various  cut-off  frequencies. 
The  scope  of  this  thesis  does  not  cover  the  multipath  effect,  so  multipath  effect 
has  been  removed  from  the  code  and  all  the  alphas  are  considered  as  one. 

Ezzat’s  solution  does  not  give  an  error  ellipsoid  for  position  estimate  error, 
nor  does  it  yield  emitter  position-covariance  as  a  function  of  TDOA  measurement 
error. 

The  simulation  of  this  thesis  is  based  on  Ezzat’s  simulation.  In  addition  to 
his  simulation,  position  estimate  error  is  based  on  multiple  simulation  runs  to  get 
position  estimate  covariance. 

The  code  of  simulation  is  modified  as  in  Appendix  B  to  decide  the  number 
of  experiments  for  minimum  estimation  error.  The  total  standard  deviation  of  the 
estimated  position  for  multiple  experiments  is  computed  in  Equation  (3.77).  The 
result  is  as  in  Figure  25. 
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Figure  25.  The  Effect  of  the  Number  of  Experiments  on  Accuracy 


It  can  be  seen  from  Figure  25  that  the  total  standard  deviation  reaches  a 
stable  level  around  500  experiments. 

The  covariance  matrix  is  based  on  the  error  resulting  from  the  noise. 

Ezzat’s  technique  does  not  show  the  effect  of  the  noise.  The  effect  of  the 
noise  to  TDOA  has  been  added  to  this  technique  and  simulation  as  in 


TDOA  =  t  +  eN 


(3.60) 


where 

t  is  the  actual  time  difference  of  arrival 
eN  is  the  error  caused  by  the  noise  (Noise  Error) 
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ewis  a  random  number  from  a  normal  distribution  with  a  0  mean  and  <jt 

standard  deviation.  Standard  deviation  is  explained  in  detail  in  Section  C  Part  2 
and  the  values  used  for  the  simulation  are  taken  from  that  part.  Because 
bandwidth,  SNR  and  integration  time  changes  for  every  receiver  pair,  noise  error 
is  different  for  each  receiver  pair,  although  the  TDOA  standard  deviation  is 
assumed  to  be  the  same  for  each  receiver  pair. 

The  simulation  is  done  for  500  times  for  efficiency  as  seen  in  Figure  25. 
Each  result  of  the  experiments  for  emitter  location  error  is  collected  in  an  error 
matrix  with  a  size  of  500*3.  500  represents  the  number  of  experiments  and  3 
represents  the  location  estimation  error  for  each  axis.  This  error  matrix  is  used  to 
compute  the  covariance  matrix  by  using  Matlab. 

Covariance  matrix  can  be  found  as 


cov(X)  =  E[(X  -  E[X])(X  -  E[X])T]  (3.61) 


If  Equation  (3.61)  is  applied  to  multiple  experiments,  the  covariance  matrix 
consists  of  variances  and  correlation  coefficients  for  every  axis  as  seen  in 
Equation  (3.62). 


'uu 


ax  Pxy°x°y  Pxz<7x(Jz 

Pxy^x^y  ®y  P yz® y^ z 

Pxz°x°z  PyzVy^z  Gz 


(3.62) 


The  covariance  matrix  shows  the  relationship  between  every  axis  (x,  y, 
and  z).  Because  the  covariance  matrix  shows  the  relationship  between  all  axis, 
we  convert  to  a  covariance  matrix  which  shows  the  only  the  variance  for  each 
single  axis  as  in 
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ax2  0  0 

0  <  0 

0  0  CTZ2 


(3.63) 


The  covariance  matrix  in  Equation  (3.63)  can  be  interpreted  as  a 
confidence  ellipsoid.  Standard  deviation  in  every  axis  of  Equation  (3.63) 
represents  half  of  the  length  of  each  axis  of  confidence  ellipsoid.  2 <jx  represents 

the  length  of  the  x-axis,  2 ay  represents  the  length  of  the  y-axis,  and  2crz 
represents  the  length  of  the  y-axis  of  the  ellipsoid. 

Because  the  covariance  matrix  in  Equation  (3.62)  shows  the  noise 
correlation  between  every  axis,  it  has  to  be  converted  to  a  diagonalized 
covariance  matrix  as  in  Equation  (3.63)  in  order  to  get  the  relationships  just  for 
axes  with  themselves,  which  can  be  interpreted  as  standard  deviations  along 
three  orthogonal  position  axes. 

Matlab  command  is  used  to  get  the  eigenvalues  of  the  covariance  matrix. 
The  eigenvalues  matrix  is  the  rotation  matrix,  which  can  be  applied  to  covariance 
matrix  in  Equation  (3.62)  to  get  the  diagonalized  covariance  matrix  in  Equation 

(3.63) .  The  following  computation  is  done  to  apply  the  rotation  matrix  to  the 
covariance  matrix  (3.62)  in  order  to  get  the  diagonalized  covariance  matrix  in 

(3.63) 


DiagonalCovMat  =  (Eigenvectors)'  *  (CovMat)'1  *  Eigenvectors  (3.64) 

The  rotation  matrix  is  explained  in  Slabaugh  and  Shoemake’s  paper 
(Slabaugh,  1999  &  Shoemake,  1985).  The  rotation  has  been  done  according  to 
the  following  order;  rotate  x  axis,  rotate  y  axis  and  rotate  z  axis. 
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R=Rz(^)Ry(#)Rx(^) 


(3.65) 


where 


Rz(^)  is  the  rotation  in  z  axis  with  the  angle  of  (f> 

Ry(#)  is  the  rotation  in  y  axis  with  the  angle  of  # 

Rx(^)  is  the  rotation  in  y  axis  with  the  angle  of  y 


Rotation  for  z  axis  is  defined  as  follows 


RZW= 


cos^ 

sin^ 


0 


-sin^  0 
cos^  0 
0  1 


Rotation  for  y  axis  is  defined  as  follows 


cos# 


Ry(*)= 


0 

sin# 


0 

1 

0 


sin# 

0 

cos# 


Rotation  for  x  axis  is  defined  as  follows 


R*(vO= 


i 

o 

0 


0 

cos  \ff 
sin^ 


0 

-siny/ 
COS  If/ 


(3.66) 


(3.67) 


(3.68) 
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As  a  result  rotation  matrix  is 


Ru  /-?12  /-?13 

^21  ^22  ^23 


R 


31 


O  P 

'  '32  '  *33 


(3.69) 


and  in  terms  of  rotation  angle 


R  = 


cos  6  cos  <j> 
costfsin^ 
-sin  0 


sin^sin0cos^-cos^sin0 
sin  y/  sin  9  sin  cj)  +  cos  y/  cos  <f> 
sin^cos^ 


cos^sin^cos^  +  sin^sin^ 
cos  y/  sin  9  sin  (f>-  sin  y/  cos  <j> 
cos  ^  cos  <9 


(3.70) 


where  the  rotation  angles  are 


<9  = -sin  ~\R3,) 


(3.71) 


f 

y/  =  atan2 

v 


p  R  \ 
r*32  r*33 

cos#’ cos 6 , 


(3.72) 


<t>  =  atan2 


»2,  i  q,  1 

vCOS#’COS#y 


(3.73) 


The  rotation  matrix  has  to  be  orthogonal;  in  other  words,  R*RT  =1,  where 

I  is  the  identity  matrix  or  R1  =RT,  where  R’1  is  inverse  of  rotation  matrix  and  RT 
is  the  transpose  of  the  rotation  matrix.  The  easiest  way  to  find  the  rotation  matrix 
is  to  find  the  eigenvectors  of  the  covariance  matrix.  A  matrix  formed  from  the 
combination  of  the  eigenvectors  of  the  covariance  matrix  is  the  rotation  matrix. 

As  explained  earlier,  the  result  in  Equation  (3.63)  gives  the  variances  for 

each  axis  diagonally.  The  square  root  of  these  variances  gives  the  standard 
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deviation  that  is  essential  to  check  the  accuracy  of  the  system.  An  ellipsoid 
represents  where  the  estimated  emitter  location  is;  the  estimated  location  of  the 
emitter  might  be  any  point  inside  that  ellipsoid.  These  variances  form  this 
ellipsoid.  Using  the  square  root  of  each  variance,  standard  deviations  represent 
the  length  of  the  half  of  each  axis.  This  ellipsoid  represents  the  volume  of  space 
which  contains  the  emitter  location  with  probability  0.25. 

A  25  percent  probability  represents  one  standard  deviation.  We  need  to 
multiply  these  standard  deviations  to  change  the  scale  of  the  ellipsoid  in  order  to 
add  different  probabilities.  Ezzat’s  algorithm  and  simulation  do  not  include  any 
probability  of  the  estimated  emitter  location. 

The  probability  is  added  to  the  simulation  as  the  length  of  the  each  semi¬ 
axis  of  the  ellipsoid  as  explained  in  Torrieri’s  paper  as  (Torrieri,  1 984) 


I x-axis  ~  4K(Jx 

=  yfc J  <3-74> 

4-„s = 


The  probability  of  the  actual  location  of  the  emitter  lies  inside  that  ellipsoid 
is  given  by  (Torrieri,  1984) 


Pe(rc)  =  erf(4K/2)-(42K/n)exp(-Kl2)  (3.75) 


where 
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(3.76) 


erf(x) 


2 


Jexp  {-t2)dt 


0 


If  the  reader  is  more  interested  in  the  statistical  theory,  he/she  can  read 
Torrieri’s  paper  (Torrieri,  1984). 


k- values  have  been  calculated  for  various  probabilities  and  shown  in  Table 


1. 


Probability 

K 

99  % 

13.89702713 

95% 

8.934259824 

90  % 

6.922544695 

85% 

5.766298681 

80  % 

4.949459443 

75% 

4.314831079 

50  % 

2.298290951 

45% 

2.008418735 

25% 

1.011546612 

Table  1 .  Probability  of  Actual  Location  of  Emitter  to  Lie  in  Ellipsoid  and  k 

Values 


The  k  values  become  larger  when  the  probability  increases,  resulting  in  a 
larger  ellipsoid  as  seen  in  Figure  26. 


Figure  26.  k  Values 
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C.  SIMULATION  RESULTS 

Since  the  standard  deviation  gives  a  good  understanding  of  the  accuracy, 
the  simulation  is  run  to  analyze  the  accuracy  with  the  standard  deviation  for 
every  axis  and  the  total  standard  deviation. 

Total  standard  deviation  is  given  by 


a total  =  ^ax2  +  Gy  +  Gz  (3-77) 

where 

o x  is  the  standard  deviation  for  x-axis 

cry  is  the  standard  deviation  for  y-axis 

o z  is  the  standard  deviation  for  z-axis 

The  simulation  gives  the  following  parameters  as  output 

•  Standard  deviation  for  each  axis  as  follows  (in  meters) 

Std_Dev_Data  = 

29.2089  (x  axis) 

30.0181  (y  axis) 

80.2068  (z  axis) 

•  Total  standard  deviation  as  follows  (in  meters) 

T  ot_Std_Dev  = 

90.4842 

•  Plot  of  the  locations  of  the  receivers  are  shown  in  Figure  27  and  an 
exaggerated  confidence  ellipsoid  is  added  in  Figure  28. 
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Figure  27.  TDOA  Geometry  in  3D 


The  exaggerated  confidence  ellipsoid  as  in  Figure  28  is  the  diagonalized 
confidence  ellipsoid  with  axes  exaggerated  by  a  factor  of  1000  for  visibility  and 
centered  on  the  estimated  position. 
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Figure  28.  TDOA  Geometry  and  Exaggerated  Confidence  Ellipsoid 

The  simulation  is  run  for  the  following  analysis 

•  The  effect  of  the  distance  on  the  accuracy, 

•  The  effect  of  the  noise  on  the  accuracy, 

•  Comparison  between  moving  receivers  and  the  stationary 
receivers, 

•  The  optimum  distribution  of  the  receivers  for  better  accuracy. 

1 .  The  Effect  of  the  Distance  on  the  Accuracy 

The  developed  simulation  is  in  Appendix  C. 

The  effect  of  the  distance  on  the  accuracy  is  analyzed  in  different  ways; 
distance  change  in  one  axis,  distance  changes  in  two  axes,  and  distance 
changes  in  all  three  axes.  Figure  29  shows  the  effect  of  the  distance  in  one  axis, 
Figure  30  shows  the  effect  of  the  distance  in  two  axes  and  Figure  31  shows  the 
effect  of  the  distances  in  three  axes. 
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Figure  29.  Effect  of  Distance  to  Accuracy  for  x-axis 
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Figure  30.  Effect  of  Distance  to  Accuracy  for  two  Axes  (x  and  y  axis) 
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Figure  31 .  Effect  of  Distance  to  Accuracy  for  Three  Axes 


The  distance  affect  the  total  standard  deviation  linearly  since  the  distance 
change  happens  in  one  or  two  axes  since  these  kind  of  changes  result  in  a 
change  in  the  shape  of  the  geometry  of  the  receiver  location.  The  distance 
change  in  two  axes  affect  the  total  standard  deviation  more  rapidly  compared  to 
the  distance  change  in  one  axis. 

On  the  other  hand,  total  standard  deviation  is  not  effected  linearly  from 
distance  change  in  three  axes  because  changes  in  three  axes  result  in  changing 
the  scale  of  the  cluster,  nothing  more.  The  average  total  standard  deviation  in 
three  axes  distance  change  is  around  300  meters  for  this  configuration  of  random 
geometry. 

2.  The  Effect  of  the  Noise  on  the  Accuracy 

Noise  affects  the  system  and  the  simulation  as 
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TDOA  -t+eN  (3.78) 

where  t  is  the  time  of  arrival  and  ewis  the  noise  error. 

The  noise  error  changes  for  each  receiver  pair.  For  this  simulation,  the 
noise  distribution  is  considered  as  a  normal  distribution  with  a  standard  deviation 
<jx  and  zero  mean  for  noise  error. 

The  developed  simulation  for  stationary  receivers  is  in  Appendix  D. 

The  simulation  is  run  for  various  standard  deviations  to  understand  the 
effect  of  the  noise  on  the  system.  Figure  32  shows  the  effect  of  the  noise  on  the 
accuracy. 


Figure  32.  The  Effect  of  Noise  on  Accuracy 


From  Figure  32,  it  can  be  seen  that  the  accuracy  in  terms  of  total  standard 
deviation  changes  linearly  with  respect  to  noise  standard  deviation.  This  is  a 
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good  method  of  increasing  the  Signal  to  Noise  Ratio  (SNR)  to  keep  the  accuracy 
high. 

This  is  consistent  with  the  theoretical  result  that  the  diagonalized 
variances  of  the  estimated  position  are  proportional  to  the  measured  TDOA 
variances. 


Figure  33.  The  Effect  of  Noise  to  Accuracy  with  Moving  UAV 

From  Figure  33,  noise  effect  with  moving  UAV  can  affect  the  accuracy  just 
like  the  way  it  affects  without  the  UAV,  but  total  standard  deviation  for  estimated 
emitter  location  with  UAV  is  smaller  than  the  one  without  the  UAV. 

3.  Comparison  Between  Moving  Receivers  to  the  Stationary  Receivers 

Ground  forces  can  use  UAVs  instead  of  stationary  receivers  even  though 
the  use  of  stationary  receivers  is  more  probable.  The  simulation  is  developed  to 
show  the  effect  of  moving  UAVs  to  the  accuracy  in  terms  of  total  standard 
deviation. 
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The  developed  simulation  is  in  Appendix  F  for  both  one  and  two  flying 

UAVs. 

The  simulation  is  run  once  for  five  random  stationary  receivers,  four 
random  stationary  receivers  plus  a  moving  UAV  and  three  random  stationary 
receivers  plus  two  moving  UAVs.  The  simulation  is  run  for  random  stationary 
receivers  first,  then  the  same  procedure  is  followed  for  four  random  stationary 
receivers  and  a  moving  UAV.  The  same  procedure  is  followed  for  the  third 
scenario  where  there  are  three  random  stationary  receivers  plus  two  moving 
UAVs.  For  each  set  of  receiver  positions,  the  emitter  is  placed  at  (0,0,0)  and  500 
TDOAs  with  random  noise  are  generated. 

The  results  are  analyzed  in  terms  of  the  standard  deviation  for  each  axis 
and  total  standard  deviation,  which  is  a  combination  of  standard  deviations  of  ail 
axes. 

Figure  34  shows  the  random  locations  of  five  stationary  receivers  which 
try  to  locate  the  emitter. 


Figure  34.  Random  Locations  of  Stationary  Receivers 

79 


Figure  35  shows  the  magnified  confidence  ellipsoid  for  five  stationary 


receivers. 


Figure  35.  Magnified  Confidence  Ellipsoid  for  Random  Stationary  Receivers 

Standard  deviations  for  each  axis  for  five  random  stationary  receivers  are: 
(tx  =  30.3517  m 
oy  =  32.7943  m 
<rz  =  77.5429  m 

Total  standard  deviation  for  five  random  stationary  receivers  is: 

(7T  =  89.4964  m 

Figure  36  shows  the  random  locations  of  four  receivers  and  a  moving  UAV 
which  try  to  locate  the  emitter. 
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The  UAV  flies  over  the  operation  area  at  a  constant  altitude  for  that 
scenario.  20  points  of  the  route  of  the  UAV  are  used  for  experiments.  500 
experiments  are  done  at  every  point  of  the  UAV.  Total  standard  deviation  is 
calculated  for  every  point  of  the  UAV.  The  results  are  collected  in  a  matrix  and 
plotted  as  in  Figure  37. 


x  10 


TDOA  Geometry 


Figure  36.  Random  Locations  of  Four  Stationary  Receivers  and  One  Moving  UAV 
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Figure  37.  Total  Standard  Deviation  Change  for  Four  Random  Stationary  Receivers 

and  a  Moving  UAV 

The  confidence  ellipsoid  for  each  experiment  is  not  plotted  but  the  numeric 
values  of  standard  deviation  for  each  axes  can  be  found  in  the  proposed 
simulation  code. 

Figure  38  shows  the  random  locations  of  three  stationary  receivers  and 
two  moving  UAV,  which  try  to  locate  the  emitter. 

The  UAVs  fly  over  the  operation  area  at  a  constant  altitude  for  that 
scenario.  20  points  of  the  routes  of  each  UAV  are  used  for  experiments.  500 
experiments  are  done  at  every  point  pair  of  the  UAVs.  Total  standard  deviation  is 
calculated  for  every  point  pair  of  the  UAVs.  The  results  are  collected  in  a  matrix 
and  plotted  as  in  Figure  39. 
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Figure  38.  Random  Locations  of  Three  Stationary  Receivers  and  Two  Moving  UAVs 


Figure  39.  Total  Standard  Deviation  Change  for  Three  Random  Stationary  Receivers 

and  two  Moving  UAVs 
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The  confidence  ellipsoid  for  each  experiment  is  not  plotted  but  the  numeric 
values  of  standard  deviation  for  each  axes  can  be  found  in  the  proposed  code. 

It  is  clear  from  these  datum  that  four  random  stationary  receivers  and  a 
flying  UAV  has  the  smaller  standard  deviation  for  each  axis,  and  they  also  have 
smaller  total  standard  deviation  compared  to  the  other  scenarios. 

4.  The  Optimum  Distribution  of  the  Receivers  for  Better  Accuracy. 

Four  receivers  can  be  positioned  around  the  emitter  using  different 
geometric  patterns,  such  as  straight,  trapezoidal,  parallelogram,  inverted  triangle, 
Y  shaped,  lozenge,  square,  and  rectangle  as  explained  in  Yan-Ping’s  article.  3D 
curves  of  position  error  are  showed  in  Figure  40  and  Figure  41 . 


Figure  40.  3D  Curves  of  Position  Error  for  Straight  (a),  Trapezoidal  (b),  Parallelogram 

(c)  and  Rectangle  (d) 
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Figure  41 .  3D  Curves  of  Position  Error  for  Lozenge  (e),  Inverted  Triangle  (f),  Square 

(g),  Y-Shaped  (h) 


Coverage  areas  are  shown  in  Table  2. 


Scene 

Coverage  value 

Straight 

10.2% 

Trapezoidal 

29.4% 

Parallelogram 

40.5% 

Rectangle 

40.9% 

Lozenge 

54.9% 

Inverted  triangle 

58.6% 

Square 

61.6% 

Y-sharp 

88.6% 

Table  2.  Coverage  Values  for  Different  Type  of  Receiver  Distributions  (From 

Yan-Ping  et  al,  2010) 
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It  is  clear  that  the  Y-shaped  receiver  distribution  is  the  best  for  coverage 
values  of  the  desired  reach  area. 

The  developed  simulation  is  in  Appendix  F. 

The  following  position  parameters  are  used  in  the  simulation  for  eight 
different  scenarios, 

Scenario  1 :  Receivers  are  around  the  emitter  (without  altitude  difference) 
as  in  Figure  42.  The  units  are  in  kilometers  (km). 
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Figure  42.  TDOA  Geometry  for  Scenario  1 
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Scenario  2:  Receivers  are  on  a  straight  line  (without  altitude  difference)  as 
in  Figure  43.  The  units  are  in  kilometers  (km). 
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Z5 

— 

Scenario  3:  Receivers  are  around  the  emitter  (with  altitude  difference)  as 
in  Figure  44  and  Figure  45.  The  units  are  in  kilometers  (km). 
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TDOA  Geometry  and  Confidence  Ellipsoid 


Figure  44.  TDOA  Geometry  for  Scenario  3  (Top  View) 


TDOA  Geometry  and  Confidence  Ellipsoid 


Figure  45.  TDOA  Geometry  for  Scenario  3  (Side  View) 
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Scenario  4:  Receivers  are  on  a  straight  line  (with  altitude  difference)  as  in 
Figure  46.  The  units  are  in  kilometers  (km). 


X 

II 

1 

M 

i 

ii 

X 

Z1  =0 

1 

II 

CM 

X 

o 

ii 

CM 

X 

Z2  =  1 

X3  =  0; 

Y3  =  2; 

Z3  =  2 

ii 

X 

Y4  =  3; 

Z4  =  3 

X5=  1.5; 

Y5  =  4; 

Z5  =  4 

TDOA  Geometry  and  Confidence  Ellipsoid 


Figure  46.  TDOA  Geometry  for  Scenario  4 

Scenario  5:  Receivers  are  in  a  trapezoidal  formation  as  in  Figure  47.  The 
units  are  in  kilometers  (km). 
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Scenario  6:  Receivers  are  in  a  parallelogram  formation  as  in  Figure  48. 
The  units  are  in  kilometers  (km). 
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Figure  48.  TDOA  Geometry  for  Scenario  6 

Scenario  7:  Receivers  are  in  a  lozenge  formation  as  in  Figure  49.  The 
units  are  in  kilometers  (km). 
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Figure  49.  TDOA  Geometry  for  Scenario  7 


Scenario  8:  Receivers  are  in  an  inverted  triangle  formation  as  in  Figure  50. 
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Table  3  shows  the  total  standard  deviations  for  all  eight  scenarios. 


Scenario 

Total  Standard 

# 

Deviation 

1 

3.9003  km 

2 

1.2128  km 

3 

67.2223  m 

4 

751.0645  m 

5 

9.5230  km 

6 

1.4135  km 

7 

1 .2795  km 

8 

3.7662  km 

Table  3.  Estimated  Emitter  Location  and  Total  Standard  Deviation  for  Position 

Scenarios 


93 


It  is  clear  from  Table  3  that  scenario  3  has  the  best  two  total  standard 
deviations.  This  type  of  receiver  positioning  should  be  used  to  get  the  best 
accuracy  if  the  system  is  using  the  proposed  close-form  geolocation  technique. 
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V.  CONCLUSION  AND  RECOMMENDATIONS 


A.  CONCLUSION 

Electronic  warfare  plays  a  dominant  role  in  today’s  technological  military 
world.  Electronic  warfare  support  is  one  of  the  core  elements  of  EW,  to  include 
direction  finding  and  geolocation. 

Passive  geolocation  of  an  emitter,  requiring  no  emission  of  energy  by  the 
receiving  and  locating  devices,  plays  an  important  role  in  both  military  and  civilian 
applications.  Passive  geolocation  can  locate  targets  by  measuring  their 
electromagnetic  emissions. 

The  greatest  internal  problem  for  the  Turkish  military  at  present  is  the 
terrorist  groups  operating  within  the  southeast  part  of  Turkey.  Military  forces  have 
been  conducting  operations  against  them  for  nearly  30  years.  The  operations 
occur  in  two  phases,  find  and  destroy.  Finding  the  location  of  the  terrorists  in  the 
mountainous  areas  of  Turkey  is  crucial.  Civilian  and  military  intelligence  has  been 
used  at  the  strategic  level  to  locate  them,  but  without  operational  and  tactical 
location  systems  success  is  fleeting  at  best. 

One  approach  that  uses  multiple  receivers  is  TDOA,  which  is  the  main 
concern  of  this  thesis.  For  this  technique,  the  time  of  arrival  of  the  received  signal 
at  different  receivers  is  measured;  the  measurements  of  the  times  of  arrival  are 
then  used  to  find  the  time  difference  of  arrival.  These  differences  are  used  to 
determine  the  location  of  the  emitter. 

Finding  the  time  difference  of  arrival  for  radar  signals  is  relatively  easy, 
compared  to  CW  signals.  Because  radar  signals  contain  pulses,  and  these 
pulses  are  almost  square  (ideally),  they  have  a  defined  beginning  and  an  end. 
These  beginnings  and  endings  provide  a  time  mark  for  time  difference  of  arrival. 
These  time  marks  are  used  for  measuring  the  times  of  arrival;  subtracting  times 
of  arrival  results  in  a  time  difference  of  arrival.  For  a  CW  signal,  there  are  no  time 
marks.  They  do  not  have  a  specific  beginning  or  end,  which  makes  time 
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difference  of  arrival  more  difficult  to  calculate.  A  cross-correlation  technique  has 
to  be  used  for  CW  signals  to  find  the  time  difference  of  arrival. 

There  is  quite  a  bit  of  literature  discussing  many  proposed  methods  to  find 
the  location  of  the  emitter  using  the  measured  time  difference  of  arrival.  Many  of 
the  techniques  require  range  data.  The  technique  discussed  in  this  thesis, 
Ezzat’s  technique,  does  not.  His  technique  is  a  closed-form  geolocation 
technique. 

The  areas  where  the  operations  take  place  are  very  mountainous. 
Because  of  this  reason,  it  is  necessary  to  use  a  3D  TDOA  technique. 

This  thesis  provided  a  brief  discussion  on  the  fundamentals  of  TDOA  and 
a  closed-form  usage  of  TDOA.  A  simulation  model  was  developed  in  order  to 
observe  the  performance  of  the  proposed  close-from  technique  in  different 
scenarios.  The  effects  of  the  distance,  noise  and  distribution  of  the  receivers  to 
accuracy  were  analyzed.  A  comparison  between  a  system  of  stationary  receivers 
and  a  system  of  combinations  of  stationary  and  moving  receivers  was  done  to 
understand  the  best  combination  of  receivers  for  best  accuracy.  For  each 
combination  of  receivers,  500  simulations  with  added  independent  noise  were 
run  to  acquire  an  estimate  of  the  emitter  position  and  to  obtain  diagonalized 
covariance  values. 

The  developed  simulation  was  based  on  Ezzat’s  closed-form  technique. 
His  technique  was  repeated  many  times  to  understand  its  robust  nature. 

Distance  affects  accuracy  in  a  linear  manner,  as  explained  in  Chapter  4 
Section  C.  Being  close  to  the  emitter  gives  the  best  accuracy.  The  closer  to  the 
emitter,  the  better  the  accuracy. 

Noise  also  has  a  linear  effect  on  accuracy,  as  explained  in  Chapter  4 
Section  C.  When  the  standard  deviation  of  the  noise  increases,  the  standard 
deviation  of  the  position  estimate  increases  proportionally,  that  results  a 
reduction  in  the  accuracy.  Since  the  noise  considered  in  the  simulation  is 
background  noise  and  cannot  be  removed  from  the  system,  the  best  way  to 
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overcome  this  problem  is  to  use  a  receiver  with  a  better  sensitivity  so  as  to 
achieve  a  higher  signal  to  noise  ratio. 

The  system  should  include  at  least  four  stationary  and  one  flying  UAV  as 
the  fifth  receiver  for  better  accuracy  as  explained  in  Chapter  4  Section  C.  It  is 
clear  from  the  simulation  that  having  a  flying  UAV  with  stationary  ground 
receivers  has  a  significant  effect  on  the  accuracy.  One  other  advantage  of  having 
a  UAV  is  not  losing  the  line  of  sight  with  the  signal  of  interest.  It  is  difficult  to  have 
line  of  sight  at  all  times,  especially  in  high  clutter  areas  that  the  Turkish  military 
forces  operate.  UAVs  can  solve  this  problem.  They  are  also  faster  and  more 
mobile  than  ground  receivers,  which  make  them  preferable  for  better  coverage. 
UAVs  also  offer  the  advantage  of  adding  the  FDOA  technique  to  the  calculation 
of  the  emitter/enemy  location,  a  discussion  that  was  beyond  the  interest  area  of 
this  thesis,  but  is  relevant  to  the  objective  of  locating  emitters. 

The  two  best  distributions  of  receivers  are  shown  in  Scenarios  3  and  4,  as 
explained  in  Chapter  4  Section  C.  One  solution  is  to  distribute  the  receivers 
around  the  emitter  with  altitude  differences  as  in  Figures  51  and  52.  Figure  51 
shows  the  geometry  from  the  side  view  and  Figure  52  shows  the  geometry  from 
the  top  view.  The  other  distribution  is  to  have  the  receivers  on  a  baseline  with 
altitude  differences  as  in  Figure  53.  The  operational  forces  should  try  to  locate 
the  receivers  as  explained  in  one  of  these  scenarios  for  best  accuracy. 
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TDOA  Geometry  and  Confidence  Ellipsoid 


Figure  51 .  TDOA  Geometry  of  the  Receivers  around  the  Emitter  (Top  View) 


TDOA  Geometry  and  Confidence  Ellipsoid 


Figure  52.  TDOA  Geometry  of  the  Receivers  around  the  Emitter  (Side  View) 
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TDOA  Geometry  and  Confidence  Ellipsoid 
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Figure  53.  TDOA  Geometry  of  the  Receivers  on  a  Straight  Line  with  Altitude 

Difference 


The  main  disadvantage  of  the  proposed  closed-form  geolocation  system  is 
a  requirement  to  use  one  more  receiver  than  with  some  other  TDOA  techniques, 
for  which  it  is  only  necessary  to  have  four  TDOA  receivers. 

The  simulation  enables  the  user  to  analyze  the  performance  of  the  closed- 
form  geolocation  technique  under  desired  conditions.  It  must  be  noted  that  the 
values  obtained  from  the  simulation  may  differ  from  real  environment  values, 
since  many  assumptions  were  made  and  certain  values  were  kept  constant  in 
order  to  simplify  the  optimization  process. 

From  the  Turkish  military  perspective,  since  the  developed  model 
examines  specific  conditions,  the  model  can  be  upgraded  to  an  advanced  level 
according  to  the  needs  and  capabilities  of  the  Turkish  military  and  can  be  an 
example  model  for  improved  future  studies. 
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B.  RECOMMENDATIONS 


EW  will  maintain  a  critical  place  on  the  battlefield.  Almost  all  nations  are 
spending  money  for  EW  research.  Every  nation  tries  to  take  part  in  the 
intelligence  arena.  All  these  facts  increase  the  importance  of  operational  tactics 
development  for  these  technologies. 

It  is  obvious  that  closed-form  geolocation  has  certain  advantages  in 
specified  areas.  The  research  conducted  under  this  thesis  tried  to  examine  these 
features  and  created  a  model  in  order  to  simulate  the  environment.  The 
simulation  model  can  be  a  starting  point  and  can  be  improved  in  several  ways  for 
future  studies. 

Contrasting  the  proposed  closed-form  geolocation  technique  with  the 
other  forms  of  TDOA  geolocation  that  use  repeated  measurements  involving  at 
least  one  moving  receiver  might  be  a  good  research  subject  to  determine  which 
TDOA  geolocation  technique  to  use  for  better  accuracy. 

Instead  of  making  assumptions,  using  fixed  values  and  omitting  certain 
facts  in  order  to  reduce  the  complexity,  allowing  for  more  detailed  studies. 
Multipath  effects  can  be  added  for  better  understanding.  Allowing  for  multipath 
may  well  produce  better  results  in  the  presence  of  mountainous  clutter.  The 
academic  studies  that  have  been  used  as  references  in  this  thesis  are  the 
primary  resources  that  can  lead  to  future  research  and  improvements  in  TDOA 
geolocation. 

Future  studies  may  try  to  answer  how  to  decrease  the  number  of  receivers 
while  maintaining  accuracy.  Because  synchronization  of  the  receivers  plays  a 
very  important  role  in  multiple  receiver  systems,  as  in  this  thesis,  future  studies 
may  focus  on  how  develop  better  synchronization  in  terms  of  system  design  and 
used  techniques. 

A  new  design  of  the  overall  system  may  be  another  area  of  research, 
where  the  TDOA  system  and  the  receivers  are  connected  to  a  central  control 
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point  over  the  Internet  or  over  a  tactical  network.  This  may  help  decision  makers 
to  develop  better  battlefield  awareness  and  more  tactical  options. 
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APPENDIX  A 


MATLAB  CODE  FOR  SIMULATION  OF  CLOSED-FORM  GEOLOCATION 

TECHNIQUE 


function  TDOA 
%  Written  by:  Volkan  TAS 

%  Based  on  code  developed  by  EZZAT  G.  BAKHOUM,  described  in 
%  Closed-Form  Solution  of  Hyperbolic  Geolocation 

%  Equations,  %  2006 

%  Latest  revision:  August  15,  2012 

%clear  the  screen 
clc 

%Variables 

%number  of  experiments  for  accuracy 
n  =  500; 

%number  of  the  positions  of  the  UAV 
m  =  1  ; 


k=8 . 934259824 ;  %k  for  95%  of  probability 


o, 

o 

k=6. 922544695; 

%k 

for 

90% 

of 

probability 

o, 

o 

k=5. 766298681; 

%k 

for 

85% 

of 

probability 

o, 

o 

k=4 . 949459443; 

%k 

for 

80% 

of 

probability 

o, 

o 

k=4 .314831079; 

%k 

for 

75% 

of 

probability 

%Coordinates  of  emitter 
E= [ 0  0  0  ]  ; 

%Speed  of  Light ( Propagation) 
c  =  3  *  1 0  A  8  ; 

%Set  All  Alphas  to  1 

alphal  =  1; 

alpha2  =  1; 

alpha3  =  1; 

alpha4  =  1; 

alpha5  =  1; 

%Stationary  Receivers 
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X3 

=-1000; 

Y3  = 

2000; 

Z3  = 

-2000; 

X4 

=  -1000; 

Y4  = 

-1000; 

Z4  = 

1000; 

X5 

=  -3000; 

Y5  = 

-1000; 

Z5  = 

-1000; 

for  i  =  1 : n, 

%initial  coordinate  matrix 
Coor2  =  [0  0  0 ] ; 

%Coordinates  of  Receivers 


for  j  = 

:  1 :  m. 

XI 

=  1000* j; 

Yl  = 

-3000* j ; 

Zl  = 

1000;  %Moving 

UAVs 

X2 

=  2000* j; 

Y2  = 

-2000* j ; 

Z2  = 

-4000* j ; 

%UAV  Coordinate  MAtrix 
XlCoor (1, j )  =  XI; 

YlCoor  (1, j )  =  Yl; 

ZlCoor  ( 1 , j )  =  Zl; 

X2Coor (1, j )  =  X2; 

Y2Coor (1, j )  =  Y2; 

Z2Coor  ( 1 , j )  =  Z2 ; 

%Calculate  the  Distance  From  emitter  to  the 
Receiver  for  TOA 

ECl=sqrt (abs (  (E (1) -XI) A2+ (E (2) -Yl) A2+  (E  (3) -Zl) A2) )  ; 
EC2=sqrt (abs (  (E (1) -X2) A2+ (E (2) -Y2) A2+ (E (3) -Z2) A2) )  ; 
EC3=sqrt (abs (  (E (1) -X3) A2+ (E (2) -Y3) A2+  (E  (3) -Z3) A2) )  ; 
EC4=sqrt (abs ( (E (1) -X4) A2+ (E (2) -Y4) A2+(E(3)-Z4)A2)); 
EC5=sqrt  (abs (  (E (1) -X5) A2+ (E  (2) -Y5) A2+ (E ( 3 ) -Z5 ) A2 ) )  ; 

%Generate  random  noise  for  every  channel  related  to 
Standat  Deviation 

Noise  Err  =  random (' Normal 0 , 5e-8 ,  1 , 5 )  ; 


%Calculate 

the  times 

of 

arrival 

tl 

=  EC1 

/ 

(alphal 

* 

c) 

+ 

Noise 

Err ( 1 ) ; 

t2 

=  EC2 

/ 

(alpha2 

* 

c) 

+ 

Noise 

~Err  (2)  ; 

t3 

=  EC3 

/ 

(alpha3 

* 

c) 

+ 

Noise 

Err ( 3 ) 

1 4 

=  EC4 

/ 

(alpha4 

* 

c) 

+ 

Noise 

Err ( 4 ) 

t5 

=  EC5 

/ 

(alpha5 

* 

c) 

+ 

Noise 

Err  (5) ; 

%from  now  on  calculate  the  emitter  location 
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%Calculate  the  coefficients  in  Eqs . (17) .  Note  that 
alpha  (assumed) , not  alpha  actual  is  used. 

all  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2/(t3-tl)  * (X3/alpha3A2  -  XI /alphal A2 ) ; 

al2  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2/(t3-tl)  *  (Y3/alpha3A2  -  Y1 /alphal A2 ) ; 

al3  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2/(t3-tl)  *  (Z3/alpha3A2  -  Zl/alphalA2) ; 

bl  =  1/ ( t2-tl )  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t3-tl)  * 

( (X3A2+Y3A2+Z3A2) /alpha3A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( t3-t2 ) ; 

a2 1  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (X4/alpha4A2  -  XI /alphal A2 ) ; 

a22  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (Y4/alpha4A2  -  Y1 /alphal A2 ) ; 

a23  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2/(t4-tl)  * (Z4/alpha4A2  -  Zl/alphalA2) ; 

b2  =  1/ ( t2-tl )  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t4-tl)  * 

( (X4A2+Y4A2+Z4A2) /alpha4A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( 1 4 - 1 2 )  ; 

a31  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2/(t5-tl)  * (X5/alpha5A2  -  XI /alphal A2 ) ; 

a32  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2/(t5-tl)  * (Y5/alpha5A2  -  Y1 /alphal A2 ) ; 

a33  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2/(t5-tl)  * (Z5/alpha5A2  -  Zl/alphalA2) ; 

b3  =  1/ (t2-tl)  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t5-tl)  * 

( (X5A2+Y5A2+Z5A2) /alpha5A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( 1 5 - 1 2 )  ; 

%Calculate  the  matrix  A  on  the  left-hand  side. 

A  =  double ([[all  al2  al3] ;  [a21  a22  a23];  [a31  a32 

a  3  3  ]  ]  )  ; 


%Calculate  the  vector  B  on  the  right-hand  side. 
B  =  double ([bl;  b2;  b3]); 


%Solve  the  simultaneous  equations  AX  =  B. 
Coorl  =  (A'*A)\A'*B; 

Coor2  =+  Coorl; 
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end 


Coor  =  Coor2/m; 

%Error  by  every  Coordinate 

Err  =  [abs (Coor  (1) -E (1) ) , abs  (Coor  (2) -E  (2) ) , 
abs  (Coor (3) — E (3) ) ] ; 

%Error  Matrix 
Mat (i, : )  =  Err; 

%Total  Error 

Total_Err  =  sqrt (Err ( 1 ) A2  +  Err(2)A2  +  Err(3)A2); 


end 

%calculate  Cov. Matrix 
Cov_Mat  =  cov (Mat); 

%Eigenvalues  and  eigenvectors 
[Eigenvectors ,  Eigenvalues  ]  =  eig (Cov_Mat ) ; 

%Find  angles  (from  "Computing  Euler  angles  from  a  rotation 
matrix" ) 

Beta  =  -asind  (Eigenvectors  ( 3 ,  1 )  )  ; 

Gamma  = 

(atan2  (Eigen_Vectors (2, 1) /cosd(Beta) , Eigen_Vectors ( 1 ,  1 )  /  cos 
d(Beta) ))* (180/pi) ; 

Alpha  = 

(atan2 (Eigen_Vectors (3,2) /cosd(Beta) , Eigen_Vectors (3, 3) /cos 
d(Beta) ))* (180/pi) ; 

%Std.  Dev  (includes  k  value) 

Std_Dev_Data  =  sqrt (diag (Eigenvalues) *k) 

%Total  Standat  Deviation 
Tot_Std_Dev  = 

sqrt (Std_Dev_Data (1) A2+Std_Dev_Data (2) A2+Std_Dev_Data (3) A2) 

%Draw  the  error  ellipsoid  from  the  diagonalized  covariance 

matrix 

hold  on 

grid  on 
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%ellipsoid 
[x,  y,  z]  = 

ellipsoid(E (1) , E (2) , E (3) , Std_Dev_Data (1) , Std_Dev_Data (2) , St 
d_Dev_Data ( 3 ) ,20) ; 

%color  of  the  ellisoid 
colormap  copper; 

%ploy  ellipsoid 
hMesh  =  mesh (x, y, z) ; 

%rotate  ellipsoid 
rotate (hMesh, [ 1  0  0], Alpha); 
rotate (hMesh, [ 0  1  0],Beta); 
rotate (hMesh, [0  0  1 ], Gamma); 

%equal  axis 
axis  equal 

%#  Change  the  camera  viewpoint 
view ([-36  18]); 

%label  axises  and  title 
xlabel ( ' X-Axis ' ) 
ylabel ( ' Y-Axis ' ) 
zlabel ( ' Z-Axis ' ) 

title ( ' TDOA  Geometry  and  Confidence  Ellipsoid') 

%plot  stationary  receivers 

plot3  (X3 , Y3 , Z3 , ' o ' ) ; 
plot3  (X4 , Y4 , Z4 , 'o' ) ; 
plot3  (X5 , Y5 , Z5 , ' o ' ) ; 

%plot  moving  receiver 

plot 3 (XlCoor, YlCoor, ZlCoor, ' ro ' ) ; 

plot3 (X2Coor, Y2Coor, Z2Coor, ' go ' ) ; 

%plot  emitter 

plot3  (E  (1) ,E  (2) ,E  (3) ,  'ms ' ) ; 
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APPENDIX  B 


MATLAB  CODE  FOR  SIMULATION  OF  CLOSED-FORM  GEOLOCATION 
TECHNIQUE  TO  DECIDE  THE  NUMBER  OF  EXPERIMENTS  FOR  OPTIMUM 

ACCURACY. 


function  TDOA_Num_Experiment_vs_Accurracy 
%  Written  by:  Volkan  TAS 

%  Based  on  code  developed  by  EZZAT  G.  BAKHOUM,  described  in 
%  Closed-Form  Solution  of  Hyperbolic  Geolocation 

%  Equations,  %  2006 

%  Latest  revision:  August  15,  2012 

%clear  the  screen 
clc 

%Variables 

%number  of  experiments  for  accuracy 
n  =  700; 

%number  of  the  positions  of  the  UAV 
m  =  1  ; 


kk=8 . 934259824; 

%  k=6. 922544695; 
%  k=5. 766298681; 
%  k=4 . 949459443; 
%  k=4 .314831079; 


%k  for  95%  of  probability 


%k 

for 

90% 

of 

probability 

%k 

for 

85% 

of 

probability 

%k 

for 

80% 

of 

probability 

%k 

for 

75% 

of 

probability 

%Coordinates  of  emitter 
E= [ 0  0  0  ]  ; 


%Speed  of  Light ( Propagation) 
c  =  3  *  1 0  A  8  ; 


%Set  All  Alphas  to  1 

alphal  =  1; 

alpha2  =  1; 

alpha3  =  1; 

alpha4  =  1; 

alpha5  =  1; 


%Stationary  Receivers 
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XI  = 

1000; 

Y1 

X2  = 

1000; 

Y2 

X3  = 

-1000; 

Y3 

X4  = 

-1000; 

Y4 

X5  = 

-3000; 

Y5 

for 

k  =  2  :  n. 

for  i  =  2 : k. 


1000; 

Z1 

=  1000; 

2000; 

Z2 

=  3000; 

2000; 

Z3 

=  -2000; 

-1000; 

Z4 

=  1000; 

-1000; 

Z5 

=  -1000; 

%Calculate  the  Distance  From  emitter  to  the 
Receiver  for  TOA 


EC1  = 

=sqrt 

(abs 

(E 

(1) 

-XI) 

A2  + 

(E 

(2) 

-Yl) 

A2  + 

(E 

(3) 

-Zl) 

A2)  )  ; 

EC2  = 

=sqrt 

(abs 

(E 

(1) 

-X2 ) 

A2  + 

(E 

(2) 

-Y2 ) 

A2  + 

(E 

(3) 

-Z2 ) 

A2)  )  ; 

EC3= 

=sqrt 

(abs 

(E 

(1) 

-X3) 

A2  + 

(E 

(2) 

-Y3) 

A2  + 

(E 

(3) 

-Z3) 

A2)  )  ; 

EC4  = 

=sqrt 

(abs 

(E 

(1) 

-X4 ) 

A2  + 

(E 

(2) 

-Y4 ) 

A2  + 

(E 

(3) 

-Z4 ) 

A2)  )  ; 

EC5= 

=sqrt 

(abs 

(E 

(1) 

-X5 ) 

A2  + 

(E 

(2) 

-Y5 ) 

A2  + 

(E 

(3) 

-Z5 ) 

A  2 )  )  ; 

%Generate  random  noise  for  every  channel  related  to 
Standat  Deviation 

Noise  Err  =  random (' Normal 0 , 5e-8 ,1,5); 


%Calculate  the  times  of  arrival. 


tl 

=  EC1 

/ 

(alphal 

* 

c) 

+ 

Noise 

Err (1) 

r 

t2 

=  EC2 

/ 

(alpha2 

* 

c) 

+ 

Noise 

Err  (2 ) 

r 

t3 

=  EC3 

/ 

(alpha3 

c) 

+ 

Noise 

"Err  (3) 

r 

t4 

=  EC4 

/ 

(alpha4 

* 

c) 

+ 

Noise 

Err (4) 

r 

t5 

=  EC5 

/ 

(alpha5 

* 

c) 

+ 

Noise 

Err  ( 5 ) 

r 

%from  now  on  calculate  the  emitter  location 

%Calculate  the  coefficients  in  Eqs . (17) .  Note  that 
alpha  (assumed) , not  alpha  actual  is  used. 

all  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2/(t3-tl)  * (X3/alpha3A2  -  XI /alphal A2 ) ; 

al2  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2/(t3-tl)  * (Y3/alpha3A2  -  Y1 /alphal A2 ) ; 

al3  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2/(t3-tl)  * (Z3/alpha3A2  -  Zl/alphalA2) ; 

bl  =  1/ (t2-tl)  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t3-tl)  * 

(  (X3A2+Y3A2  +  Z3A2) /alpha3A2-  (X1A2+Y1A2  +  Z1A2) /alphalA2)  + 
cA2  *  ( t3-t2 ) ; 

a2 1  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (X4/alpha4A2  -  XI /alphal A2); 


110 


a22  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (Y4/alpha4A2  -  Y1 /alphal A2 ) ; 

a23  =  2/(t2-tl)  *  ( Z 2 / alpha2 A2  -  Zl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (Z4/alpha4A2  -  Zl/alphalA2) ; 

b2  =  1/ ( t2-tl )  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t4-tl)  * 

( (X4A2+Y4A2+Z4A2) /alpha4A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  (t4-t2); 

a31  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2/(t5-tl)  * (X5/alpha5A2  -  XI /alphal A2 ) ; 

a32  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2/(t5-tl)  * (Y5/alpha5A2  -  Y1 /alphal A2 ) ; 

a33  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2/(t5-tl)  * (Z5/alpha5A2  -  Zl/alphalA2) ; 

b3  =  1/ (t2-tl)  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t5-tl)  * 

( (X5A2+Y5A2+Z5A2) /alpha5A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( t5-t2 ) ; 

%Calculate  the  matrix  A  on  the  left-hand  side. 

A  =  double ([[all  al2  al3] ;  [a21  a22  a23];  [a31  a32 

a  3  3  ]  ]  )  ; 

%Calculate  the  vector  B  on  the  right-hand  side. 

B  =  double ([bl;  b2;  b3]); 

%Solve  the  simultaneous  equations  AX  =  B. 

Coor  =  (A' *A) \A ' *B; 

%Error  by  every  Coordinate 

Err  =  [abs (Coor  ( 1 ) -E  ( 1 ) ) , abs  (Coor  (2 ) -E  (2 ) ) , 
abs  (Coor  (3) -E  (3) ) ] ; 

%Error  Matrix 
Mat  (i,  : )  =  Err; 


end 

%calculate  Cov. Matrix 
Cov_Mat  =  cov (Mat); 

%Eigenvalues  and  eigenvectors 
[Eigenvectors ,  Eigenvalues  ]  =  eig  (Cov_Mat )  ; 

%Std.  Dev  (includes  k  value) 

Std_Dev_Data  =  sqrt (diag (Eigenvalues) *kk)  ; 
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%Total  Standat  Deviation 
Tot_Std_Dev  = 

sqrt (Std_Dev_Data (1) A2+Std_Dev_Data (2) A2+Std_Dev_Data (3) A2) 


Tot  Std  Devi (k, :)  =  Tot  Std  Dev; 


end 

N  =  linspace  ( 1 , n, n) ; 

%plot  the  result 
plot  (N, Tot_Std_Devl ) 

%label  axises  and  title 

xlabel (' Number  of  Experiments') 

ylabel  (' Total  Standard  Deviation  (m)  ') 

title  (' Effect  of  Number  of  Experiment  to  the  Accuracy') 
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APPENDIX  C 


MATLAB  CODE  FOR  SIMULATION  OF  CLOSED-FORM  GEOLOCATION 
TECHNIQUE  TO  SEE  THE  EFFECTS  OF  DISTANCE  CHANGES  ON 

ACCURACY. 

function  TDOA_Distance_to_Accuracy 
%  Written  by:  Volkan  TAS 

%  Based  on  code  developed  by  EZZAT  G.  BAKHOUM,  described  in 
%  Closed-Form  Solution  of  Hyperbolic  Geolocation 
%  Equations,  %  2006 
%  Latest  revision:  August  15,  2012 

%clear  the  screen 
clc 

%Variables 


%number  of  experiments  for  accuracy 
n  =  500; 


kk=8 . 934259824 ;  %k  for  95%  of  probability 


k=6. 922544695;  %k  for  90% 
k=5. 766298681;  %k  for  85% 
k=4 . 949459443;  %k  for  80% 
k=4 . 314831079;  %k  for  75% 


of  probability 
of  probability 
of  probability 
of  probability 


%Coordinates  of  emitter 
E= [ 0  0  0  ]  ; 

%Speed  of  Light ( Propagation) 
c  =  3  *  1 0  A  8  ; 

%Set  All  Alphas  to  1 

alphal  =  1; 

alpha2  =  1; 

alpha3  =  1; 

alpha4  =  1; 

alpha5  =  1; 

%distance  multipication  constant 
k= 1 0  0 ; 


%Distance  Change  loop 
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for  1 


1 :  k 


%Distance  Change  in  three  axis 


%Stationary  Rece 

ivers 

XI 

=  1000*1; 

Y1  = 

1000*1; 

Z1  = 

=  1000*1; 

X2 

=  1000*1; 

Y2  = 

2000*1; 

Z2  = 

=  3000*1; 

X3 

=-1000*1; 

Y3  = 

2000*1; 

Z3  = 

=  -2000*1 

X4 

=  -1000*1; 

Y4  = 

-1000*1; 

Z4  = 

=  1000*1; 

X5 

=  -3000*1; 

Y5  = 

-1000*1; 

Z5  = 

=  -1000*1 

o, 

O 

Distance  Change  in  two 

axis 

O, 

O 

%Stationary  Receivers 

Q, 

O 

XI  =  1000; 

Y1  = 

1000*1; 

Z1 

=  1000*1; 

O, 

O 

X2  =  1000; 

Y2  = 

2000*1; 

Z2 

=  3000*1; 

O, 

O 

X3  =-1000; 

Y3  = 

2000*1; 

Z3 

=  -2000*1; 

o, 

o 

X4  =  -1000; 

Y4  = 

-1000*1; 

Z4 

=  1000*1; 

Q, 

O 

O, 

X5  =  -3000; 

Y5  = 

-1000*1; 

Z5 

=  -1000*1; 

O 

O, 

O 

Distance  Change  in  one 

axis 

O, 

O 

Stationary  Receivers 

o, 

o 

XI  =  1000; 

Y1  = 

1000; 

Z1  = 

1000*1; 

o, 

o 

X2  =  1000; 

Y2  = 

2000; 

Z2  = 

3000*1; 

Q. 

O 

X3  =-1000; 

Y3  = 

2000; 

Z3  = 

-2000*1; 

Q. 

O 

X4  =  -1000; 

Y4  = 

-1000; 

Z4  = 

1000*1; 

O 

O 

X5  =  -3000; 

Y5  = 

-1000; 

Z5  = 

-1000*1; 

for  i  =  2 : n, 

%Calculate  the  Distance  From  emitter  to  the 
Receiver  for  TOA 


EC1 

=sqrt 

(abs  ( 

(E  (1 

)  -xi) 

A2  + 

(E 

(2) 

-  Y1 

)  A2  + 

(E 

(3 

)  -zi) 

A2 ) 

)  ; 

EC2 

=sqrt 

(abs  ( 

(E  (1 

)  -X2) 

A2  + 

(E 

(2) 

-Y2 

)  A2  + 

(E 

(3 

)  -Z2) 

A2 ) 

)  ; 

EC3 

=sqrt 

(abs  ( 

(E  (1 

)  -X3) 

A2  + 

(E 

(2) 

-Y3 

)  A2  + 

(E 

(3 

)  —  Z  3 ) 

A2 ) 

)  ; 

EC4 

=sqrt 

(abs  ( 

(E  (1 

)  -X4) 

A2  + 

(E 

(2) 

-Y4 

)  A2  + 

(E 

(3 

)  —  Z  4  ) 

A2 ) 

)  ; 

EC5 

=sqrt 

(abs  ( 

(E  (1 

)  -X5) 

A2  + 

(E 

(2) 

-Y5 

)  A2  + 

(E 

(3 

)  —  Z  5 ) 

A2 ) 

)  ; 

%Generate  random 

noise 

fo 

r 

eve 

ry 

channel 

related 

to 

Standat  Deviation 

Noise  Err  =  random (' Normal 0 , 5e-8 ,1,5); 


%Calculate 

the  times 

of 

arrival 

. 

tl 

=  EC1 

/ 

(alphal 

* 

c) 

+ 

Noise 

Err 

(1) 

r 

t2 

=  EC2 

/ 

(alpha2 

* 

c) 

+ 

Noise 

Err 

(2) 

r 

t3 

=  EC3 

/ 

(alpha3 

* 

c) 

+ 

Noise 

Err 

(3) 

r 

1 4 

=  EC4 

/ 

(alpha4 

* 

c) 

+ 

Noise 

Err 

(4) 

r 
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t5  =  EC5  /  (alpha5  *  c)  +  Noise_Err ( 5 ) ; 


%from  now  on  calculate  the  emitter  location 


%Calculate  the  coefficients  in  Eqs . (17) .  Note  that 
alpha  (assumed) , not  alpha  actual  is  used. 

all  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2/(t3-tl)  * (X3/alpha3A2  -  XI /alphal A2 ) ; 

al2  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2/(t3-tl)  * (Y3/alpha3A2  -  Y1 /alphal A2 ) ; 

al3  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2/(t3-tl)  * (Z3/alpha3A2  -  Zl/alphalA2) ; 

bl  =  1/ (t2-tl)  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t3-tl)  * 

( (X3A2+Y3A2+Z3A2) /alpha3A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( t3-t2 ) ; 

a2 1  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (X4/alpha4A2  -  XI /alphal A2 ) ; 

a22  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (Y4/alpha4A2  -  Y1 /alphal A2 ) ; 

a23  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (Z4/alpha4A2  -  Zl/alphalA2) ; 

b2  =  1/ (t2-tl)  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t4-tl)  * 

( (X4A2+Y4A2+Z4A2) /alpha4A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( 1 4 - 1 2 )  ; 

a31  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2/(t5-tl)  * (X5/alpha5A2  -  XI /alphal A2 ) ; 

a32  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2/(t5-tl)  * (Y5/alpha5A2  -  Y1 /alphal A2 ) ; 

a33  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2/(t5-tl)  * (Z5/alpha5A2  -  Zl/alphalA2) ; 

b3  =  1/  ( t2-tl )  *  ( (X2A2+Y2A2  +  Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t5-tl)  * 

(  (X5A2+Y5A2  +  Z5A2) /alpha5A2-  (X1A2+Y1A2  +  Z1A2) /alphalA2)  + 
cA2  *  ( t5-t2 )  ; 

%Calculate  the  matrix  A  on  the  left-hand  side. 

A  =  double ([[all  al2  al3] ;  [a21  a22  a23];  [a31  a32 

a  3  3  ]  ]  )  ; 


%Calculate  the  vector  B  on  the  right-hand  side. 
B  =  double ([bl;  b2;  b3]); 
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%Solve  the  simultaneous  equations  AX  =  B. 
Coor  =  (A' *A) \A ' *B; 

%Error  by  every  Coordinate 

Err  =  [abs (Coor  ( 1 ) -E  ( 1 ) ) , abs  (Coor  (2 ) -E  (2 ) ) , 
abs  (Coor  (3) -E  (3) ) ] ; 

%Error  Matrix 

Mat  (i,  : )  =  Err; 


end 

%calculate  Cov. Matrix 
Cov_Mat  =  cov (Mat); 

%Eigenvalues  and  eigenvectors 
[Eigenvectors ,  Eigenvalues  ]  =  eig  (Cov_Mat )  ; 

%Std.  Dev  (includes  k  value) 

Std_Dev_Data  =  sqrt (diag (Eigenvalues) *kk) ; 

%Total  Standat  Deviation 
Tot_Std_Dev  = 

sqrt (Std_Dev_Data (1) A2+Std_Dev_Data (2) A2+Std_Dev_Data (3) A2) 


Tot  Std  Devl(l,:)  =  Tot  Std  Dev; 


end 

N  =  linspace ( 1 , k, k) ; 

%plot  the  result 
plot  (N, Tot_Std_Devl ) 

%label  axises  and  title 

xlabel (' Distance  Multiplier') 

ylabel  (' Total  Standard  Deviation  (m)  ') 

title  (' Effect  of  Distance  to  the  Accuracy') 
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APPENDIX  D 


MATLAB  CODE  FOR  SIMULATION  OF  CLOSED-FORM  GEOLOCATION 
TECHNIQUE  WITH  STATIONARY  RECEIVERS  TO  SEE  THE  EFFECTS  OF 

NOISE  ON  ACCURACY. 

function  TDOA_Noise_vs_Accuracy_Stationary_Receivers 
%  Written  by:  Volkan  TAS 

%  Based  on  code  developed  by  EZZAT  G.  BAKHOUM,  described  in 
%  Closed-Form  Solution  of  Hyperbolic  Geolocation 

%  Equations,  %  2006 

%  Latest  revision:  August  15,  2012 

%clear  the  screen 
clc 

%Variables 

%number  of  experiments  for  accuracy 
n  =  500; 

%standard  deviation  of  noise  chage  constant 
1=2  5; 


%  k  values  for  probability 

kk=8 . 934259824 ;  %k  for  95%  of  probability 


o, 

o 

kk=6. 922544695; 

%k 

for 

90% 

of 

probability 

o 

o 

kk=5. 766298681; 

%k 

for 

85% 

of 

probability 

g, 

o 

kk=4 . 949459443; 

%k 

for 

80% 

of 

probability 

g, 

o 

kk=4 .314831079; 

%k 

for 

75% 

of 

probability 

%Coordinates  of  emitter 
E= [ 0  0  0  ]  ; 

%Speed  of  Light ( Propagation) 
c  =  3  *  1 0  A  8 ; 

%Set  All  Alphas  to  1 

alphal  =  1; 

alpha2  =  1; 

alpha3  =  1; 

alpha4  =  1; 

alpha5  =  1; 
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%Stationary 

Receivers 

XI 

=  1000; 

Yl 

=  1000; 

Zl  = 

=  1000; 

X2 

=  1000; 

Y2 

=  2000; 

Z2  = 

=  3000; 

X3 

=-1000; 

Y3 

=  2000; 

Z3  = 

=  -2000; 

X4 

=  -1000; 

Y4 

=  -1000; 

Z4  = 

=  1000; 

X5 

=  -3000; 

Y5 

=  -1000; 

Z5  = 

=  -1000; 

%for  loop  for  changing  the  standard  deviation 
for  k  =  1:1, 

%for  loop  for  experiment  for  every  standard  deviation 
for  i  =  1 : n, 

%Calculate  the  Distance  From  emitter  to  the 
Receiver  for  TOA 

ECl=sqrt (abs (  (E (1) -XI) A2+ (E (2) -Yl) A2+ (E (3) -Zl) A2) )  ; 
EC2=sqrt (abs (  (E (1) -X2) A2+ (E  (2) -Y2) A2+  (E  ( 3 ) -Z2 ) A2 ) )  ; 
EC3=sqrt (abs ( (E (1) -X3) A2+ (E (2) -Y3) A2+ (E (3) -Z3) A2) ) ; 
EC4=sqrt (abs ( (E (1) -X4) A2+ (E (2) -Y4) A2+ (E ( 3 ) -Z4 ) A2 ) ) ; 
EC5=sqrt  (abs  (  (E  ( 1 )  -X5 )  A2+  (E  (2)  -Y5)  A2+(E(3)-Z5)A2)); 

%Generate  random  noise  for  every  channel  related  to 
Standat  Deviation 

Noise  Err  =  random (' Normal 0 , k* le-8 ,  1 , 5 )  ; 


%Calculate 

the  times 

of 

arrival 

• 

tl 

=  EC1 

/ 

(alphal 

* 

c) 

+ 

Noise 

_Err  (1) 

t2 

=  EC2 

/ 

(alpha2 

* 

c) 

+ 

Noise 

Err  (2 )  ; 

t3 

=  EC3 

/ 

(alpha3 

-k 

c) 

+ 

Noise 

“Err  (3) 

t4 

=  EC4 

/ 

(alpha4 

* 

c) 

+ 

Noise 

Err (4) 

t5 

=  EC5 

/ 

(alpha5 

* 

c) 

+ 

Noise 

Err  (5) ; 

%from  now  on  calculate  the  emitter  location 

%Calculate  the  coefficients  in  Eqs . (17) .  Note  that 
alpha  (assumed) , not  alpha  actual  is  used. 

all  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2/(t3-tl)  * (X3/alpha3A2  -  XI /alphal A2 ) ; 

al2  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2/(t3-tl)  *  (Y3/alpha3A2  -  Yl /alphal A2 ) ; 

al3  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2/(t3-tl)  * (Z3/alpha3A2  -  Zl/alphalA2) ; 

bl  =  1/ (t2-tl)  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t3-tl)  * 
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( (X3A2+Y3A2+Z3A2) /alpha3A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( 1 3 - 1 2 ) ; 

a2 1  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (X4/alpha4A2  -  XI /alphal A2 ) ; 

a22  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (Y4/alpha4A2  -  Y1 /alphal A2 ) ; 

a23  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (Z4/alpha4A2  -  Zl/alphalA2) ; 

b2  =  1/ (t2-tl)  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t4-tl)  * 

( (X4A2+Y4A2+Z4A2) /alpha4A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( 1 4 - 1 2 )  ; 

a31  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2/(t5-tl)  * (X5/alpha5A2  -  XI /alphal A2 ) ; 

a32  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2/(t5-tl)  * (Y5/alpha5A2  -  Y1 /alphal A2 ) ; 

a33  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2/(t5-tl)  * (Z5/alpha5A2  -  Zl/alphalA2) ; 

b3  =  1/ (t2-tl)  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t5-tl)  * 

(  (X5A2+Y5A2  +  Z5A2) /alpha5A2-  (X1A2+Y1A2  +  Z1A2) /alphalA2)  + 
cA2  *  ( t5-t2 )  ; 

%Calculate  the  matrix  A  on  the  left-hand  side. 

A  =  double ([[all  al2  al3] ;  [a21  a22  a23];  [a31  a32 

a  3  3  ]  ]  )  ; 


%Calculate  the  vector  B  on  the  right-hand  side. 
B  =  double ([bl;  b2;  b3]); 

%Solve  the  simultaneous  equations  AX  =  B. 

Coor  =  (A' *A) \A ' *B; 

%Error  by  every  Coordinate 

Err  =  [abs (Coor (1) -E (1) ) , abs  (Coor  (2) -E  (2) ) , 
abs  (Coor  (3) -E  (3) ) ] ; 

%Error  Matrix 
Mat (i, : )  =  Err; 

end  %i 

%calculate  Cov. Matrix 
Cov  Mat  =  cov (Mat); 
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%Eigenvalues  and  eigenvectors 
[Eigenvectors ,  Eigenvalues ]  =  eig  (Cov_Mat )  ; 

%Std.  Dev  (includes  k  value) 

Std_Dev_Data  =  sqrt (diag (Eigenvalues) *kk) ; 

%Total  Standat  Deviation 
Tot_Std_Dev  = 

sqrt (Std_Dev_Data (1) A2+Std_Dev_Data (2) A2+Std_Dev_Data (3)  A2) 


Tot  Std  Devi (k, :)  =  Tot  Std  Dev; 


end  %k 


N  =  linspace  (1,1,1); 


%plot  the  result 
plot  (N, Tot_Std_Devl ) 


%label  axises 
xlabel ( 'Noise 
ylabel ( ' Total 
title  (  ' Effect 


and  title 
Standard 
Standard 
of  Noise 


Deviation  (sec) ') 
Deviation  (m) ' ) 
to  the  Accuracy') 


120 


APPENDIX  E 


MATLAB  CODE  FOR  SIMULATION  OF  CLOSED-FORM  GEOLOCATION 
TECHNIQUE  WITH  ONE  MOVING  UAV  TO  SEE  THE  EFFECTS  OF  NOISE  ON 

ACCURACY. 

function  TDOA_Noise_vs_Accuracy_With_Flying_UAV 
%  Written  by:  Volkan  TAS 

%  Based  on  code  developed  by  EZZAT  G.  BAKHOUM,  described  in 
%  Closed-Form  Solution  of  Hyperbolic  Geolocation 

%  Equations,  %  2006 

%  Latest  revision:  August  15,  2012 

%clear  the  screen 
clc 

%Variables 

%number  of  experiments  for  accuracy 
n  =  500; 

Mat  =  zeros (n, 3 ) ; 

%number  of  the  positions  of  the  UAV 
m  =  10; 

%standard  deviation  of  noise  chage  constant 
1=2  5; 


%  k  values  for  probability 

kk=8 . 934259824 ;  %k  for  95%  of  probability 


o, 

o 

kk=6. 922544695; 

%k 

for 

90% 

of 

probability 

o, 

o 

kk=5. 766298681; 

%k 

for 

85% 

of 

probability 

o 

o 

kk=4 . 949459443; 

%k 

for 

80% 

of 

probability 

o, 

o 

kk=4 .314831079; 

%k 

for 

75% 

of 

probability 

%Coordinates  of  emitter 
E= [ 0  0  0  ]  ; 

%Speed  of  Light ( Propagation) 
c  =  3  *  1 0  A  8  ; 

%Set  All  Alphas  to  1 
alphal  =  1; 
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alpha2  =  1; 
alpha3  =  1; 
alpha4  =  1; 
alpha5  =  1; 


%Stationary 

Receivers 

X2  =  1000; 

Y2 

=  2000; 

Z2  = 

3000; 

X3  =-1000; 

Y3 

=  2000; 

Z3  = 

-2000; 

X4  =  -1000; 

Y4 

=  -1000; 

Z4  = 

1000; 

X5  =  -3000; 

Y5 

=  -1000; 

Z5  = 

-1000; 

%for  loop  for  changing  the  standard  deviation 
for  k  =  1:1, 


%for  loop  for  experiment  for  every  standard  deviation 
for  i  =  1 : n, 

%initial  coordinate  matrix 
Coor2  =  [0  0  0 ] ; 

%for  loop  for  the  positions  of  the  UAV 
for  j  =  l:m, 

XI  =  1000*j;  Y1  =  1000*j;  Z1  =  1000; 

%Moving  UAV 


Receiver 

Zl) A2 ) ) ; 
Z  2  )  A  2  )  )  ; 
Z3) A2 ) ) ; 
Z4  )  A2 ) )  ; 
Z5 ) A2 ) )  ; 


%Calculate  the  Distance  From  emitter  to  the 
for  TOA 


EC1  = 

=sqrt 

(abs  ( 

(E 

(1) 

-XI) 

A2  + 

(E 

(2) 

-Yl) 

A2  + 

(E 

(3)- 

EC2  = 

=sqrt 

(abs  ( 

(E 

(1) 

-X2) 

A2  + 

(E 

(2) 

-Y2 ) 

A2  + 

(E 

(3)- 

EC3= 

=sqrt 

(abs  ( 

(E 

(1) 

-X3) 

A2  + 

(E 

(2) 

-Y3) 

A2  + 

(E 

(3)- 

EC4  = 

=sqrt 

(abs  ( 

(E 

(1) 

-X4 ) 

A2  + 

(E 

(2) 

-Y4 ) 

A2  + 

(E 

(3)- 

EC5= 

=sqrt 

(abs  ( 

(E 

(1) 

-X5 ) 

A2  + 

(E 

(2) 

-Y5 ) 

A2  + 

(E 

(3)- 

%Generate  random  noise  for  every  channel 
related  to  Standat  Deviation 

Noise_Err  =  random (' Normal 0 , k* le-8 , 1 , 5 ) ; 

%Calculate  the  times  of  arrival. 

tl  =  EC1  /  (alphal  *  c)  +  Noise_Err ( 1 ) ; 

t2  =  EC2  /  (alpha2  *  c)  +  Noise_Err (2 ) ; 
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t3  =  EC3  /  (alpha3  *  c)  +  Noise_Err ( 3 ) ; 
t4  =  EC4  /  (alpha4  *  c)  +  Noise_Err ( 4 ) ; 
t5  =  EC5  /  (alpha5  *  c)  +  Noise_Err ( 5 ) ; 


%from  now  on  calculate  the  emitter  location 


%Calculate  the  coefficients  in  Eqs . (17) .  Note 
that  alpha  (assumed) , not  alpha  actual  is  used. 

all  =  2 /  ( t2-t 1 )  *  (X2/alpha2A2  -  Xl/alphalA2) 
2/(t3-tl)  * (X3/alpha3A2  -  XI /alphal A2 ) ; 

al2  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2) 
2/(t3-tl)  *  (Y3/alpha3A2  -  Y1 /alphal A2 ) ; 

al3  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2) 
2/(t3-tl)  *  (Z3/alpha3A2  -  Zl/alphalA2) ; 

bl  =  1/ ( t2-tl )  *  (  (X2A2+Y2A2  +  Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t3-tl)  * 

( (X3A2+Y3A2+Z3A2) /alpha3A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( 1 3 - 1 2 ) ; 

a21  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2) 

2 /  ( 1 4-t 1 )  * (X4/alpha4A2  -  XI /alphal A2 ) ; 

a22  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2) 

2 /  ( 1 4-t 1 )  * (Y4/alpha4A2  -  Y1 /alphal A2 ) ; 

a23  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2) 

2 /  ( 1 4-t 1 )  * (Z4/alpha4A2  -  Zl/alphalA2) ; 

b2  =  1/ ( t2-tl )  *  (  (X2A2+Y2A2  +  Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t4-tl)  * 

( (X4A2+Y4A2+Z4A2) /alpha4A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( 1 4 - 1 2 )  ; 

a31  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2) 
2/(t5-tl)  * (X5/alpha5A2  -  XI /alphal A2 ) ; 

a32  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2) 
2/(t5-tl)  * (Y5/alpha5A2  -  Y1 /alphal A2 ) ; 

a33  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2) 
2/(t5-tl)  * (Z5/alpha5A2  -  Zl/alphalA2) ; 

b3  =  1/ (t2-tl)  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t5-tl)  * 

( (X5A2+Y5A2+Z5A2) /alpha5A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( t5-t2 ) ; 

%Calculate  the  matrix  A  on  the  left-hand  side. 
A  =  double ([[all  al2  al3] ;  [a21  a22  a23];  [a31 

a32  a33 ] ] ) ; 


%Calculate  the  vector  B  on  the  right-hand  side 
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B  =  double ( [bl;  b2;  b3]); 

%Solve  the  simultaneous  equations  AX  =  B. 
Coorl  =  (A' *A) \A ' *B; 

Coor2  =+  Coorl; 
end  % j 

Coor  =  Coor2/m; 

%Error  by  every  Coordinate 

Err  =  [abs (Coor (1) -E (1) ) , abs  (Coor  (2) -E  (2) ) , 
abs  (Coor  (3) -E  (3) ) ] ; 

%Error  Matrix 

Mat  (i,  : )  =  Err; 


end  %i 

%calculate  Cov. Matrix 
Cov_Mat  =  cov (Mat); 

%Eigenvalues  and  eigenvectors 
[Eigenvectors ,  Eigenvalues  ]  =  eig  (Cov_M&t )  ; 


%Std.  Dev  (includes  k  value) 

Std_Dev_Data  =  sqrt (diag (Eigenvalues) *kk) ; 

%Total  Standat  Deviation 
Tot_Std_Dev  = 

sqrt (Std_Dev_Data (1) A2+Std_Dev_Data (2) A2+Std_Dev_Data (3) A2) 


Tot_Std_Devl ( k, : )  =  Tot_Std_Dev; 
end  %k 

N  =  linspace  ( 1 , 1, 1) ; 

%plot  the  result 
plot  (N, Tot_Std_Devl ) 
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%label  axises 
xlabel ( 'Noise 
ylabel ( ' Total 
title  (  ' Effect 


and  title 
Standard 
Standard 
of  Noise 


Deviation  (sec))') 
Deviation  (m) ' ) 
to  the  Accuracy') 
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APPENDIX  F 


MATLAB  CODE  FOR  SIMULATION  OF  CLOSED-FORM  GEOLOCATION 
TECHNIQUE  TO  SEE  THE  EFFECTS  OF  HAVING  A  MOVING  UAV  ALONG 
WITH  STAIONARY  RECEIVERS  ON  ACCURACY. 

function  TDOA_Flying_UAV 
%  Written  by:  Volkan  TAS 

%  Based  on  code  developed  by  EZZAT  G.  BAKHOUM,  described  in 
%  Closed-Form  Solution  of  Hyperbolic  Geolocation 
%  Equations,  %  2006 
%  Latest  revision:  August  26,  2012 

%clear  the  screen 
clc 

%Variables 

%number  of  experiments  for  accuracy 
n  =  500; 

Mat  =  zeros (n, 3 ) ; 

%number  of  the  positions  of  the  UAV 
m  =  2  0; 


k=8 . 934259824 ;  %k  for  95%  of  probability 


o, 

o 

k=6. 922544695; 

%k 

for 

90% 

of 

probability 

o 

o 

k=5. 766298681; 

%k 

for 

85% 

of 

probability 

o, 

o 

k=4 . 949459443; 

%k 

for 

80% 

of 

probability 

o, 

o 

k=4 .314831079; 

%k 

for 

75% 

of 

probability 

%Coordinates  of  emitter 
E= [ 0  0  0  ]  ; 

%Speed  of  Light ( Propagation) 
c  =  3  *  1 0  A  8  ; 

%Set  All  Alphas  to  1 

alphal  =  1; 

alpha2  =  1; 

alpha3  =  1; 

alpha4  =  1; 

alpha5  =  1; 
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%Stationary  Receivers 


X2 

=  1000; 

Y2  = 

-2000; 

Z2  = 

-2100; 

X3 

=-1000; 

Y3  = 

2000; 

Z3  = 

-2000; 

X4 

=  -1000; 

Y4  = 

-1000; 

Z4  = 

1000; 

X5 

=  -3000; 

Y5  = 

-1000; 

Z5  = 

-1000; 

%for  loop  for  UAV  movement 
for  j  =  l:m, 

%locations  of  the  UAV 

XI  =  1 0  0  0  *  j ;  Y1  =  1 5  0  0  *  j ;  Z1  =  2000;  %Moving  UAV 


%UAV  Coordinate  MAtrix 
XlCoor (1, j )  =  XI; 

YlCoor  (1, j )  =  Yl; 

ZlCoor ( 1 , j )  =  Zl; 

%for  loop  for  experiments  for  every  location  of  the  UAV 
for  i  =  1 ; n, 

%initial  coordinate  matrix 
Coor2  =  [0  0  0 ] ; 

%Calculate  the  Distance  From  emitter  to  the 
Receiver  for  TOA 

ECl=sqrt (abs (  (E  (1) -XI) A2+ (E (2) -Yl) A2+ (E ( 3 ) -Z 1 ) A2 ) )  ; 
EC2=sqrt (abs ( (E (1) -X2) A2+ (E (2) -Y2) A2+ (E ( 3 ) -Z2 ) A2 ) ) ; 
EC3=sqrt (abs ( (E (1) -X3) A2+ (E (2) -Y3) A2+ (E (3) -Z3) A2) ) ; 
EC4=sqrt (abs (  (E (1) -X4) A2+ (E (2) -Y4) A2+ (E ( 3 ) -Z4 ) A2 ) )  ; 
EC5=sqrt (abs ( (E (1) -X5) A2+ (E (2) -Y5) A2+ (E (3) -Z5) A2) ) ; 

%Generate  random  noise  for  every  channel  related  to 
Standat  Deviation 

Noise  Err  =  random (' Normal 0 , 5e-8 , 1 , 5 ) ; 


%Calculate  the  times  of  arrival. 


tl 

=  EC1 

/ 

(alphal 

* 

c) 

+ 

Noise 

Err (1) 

r 

t2 

=  EC2 

/ 

(alpha2 

* 

c) 

+ 

Noise 

Err  (2 ; 

r 

t3 

=  EC3 

/ 

(alpha3 

* 

c) 

+ 

Noise 

Err ( 3 ] 

r 

t4 

=  EC4 

/ 

(alpha4 

* 

c) 

+ 

Noise 

Err (4) 

r 

t5 

=  EC5 

/ 

(alpha5 

* 

c) 

+ 

Noise 

Err  ( 5 ) 

r 

%from  now  on  calculate  the  emitter  location 
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%Calculate  the  coefficients  in  Eqs . (17) .  Note  that 
alpha  (assumed) , not  alpha  actual  is  used. 

all  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2/(t3-tl)  * (X3/alpha3A2  -  XI /alphal A2 ) ; 

al2  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2/(t3-tl)  *  (Y3/alpha3A2  -  Y1 /alphal A2 ) ; 

al3  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2/(t3-tl)  *  (Z3/alpha3A2  -  Zl/alphalA2) ; 

bl  =  1/ ( t2-tl )  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t3-tl)  * 

( (X3A2+Y3A2+Z3A2) /alpha3A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( t3-t2 ) ; 

a2 1  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (X4/alpha4A2  -  XI /alphal A2 ) ; 

a22  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (Y4/alpha4A2  -  Y1 /alphal A2 ) ; 

a23  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2/(t4-tl)  * (Z4/alpha4A2  -  Zl/alphalA2) ; 

b2  =  1/ ( t2-tl )  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t4-tl)  * 

( (X4A2+Y4A2+Z4A2) /alpha4A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( 1 4 - 1 2 )  ; 

a31  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2/(t5-tl)  * (X5/alpha5A2  -  XI /alphal A2 ) ; 

a32  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2/(t5-tl)  * (Y5/alpha5A2  -  Y1 /alphal A2 ) ; 

a33  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2/(t5-tl)  * (Z5/alpha5A2  -  Zl/alphalA2) ; 

b3  =  1/ (t2-tl)  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t5-tl)  * 

( (X5A2+Y5A2+Z5A2) /alpha5A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( 1 5 - 1 2 )  ; 

%Calculate  the  matrix  A  on  the  left-hand  side. 

A  =  double ([[all  al2  al3] ;  [a21  a22  a23];  [a31  a32 

a  3  3  ]  ]  )  ; 


%Calculate  the  vector  B  on  the  right-hand  side. 
B  =  double ([bl;  b2;  b3]); 


%Solve  the  simultaneous  equations  AX  =  B. 
Coorl  =  (A'*A)\A'*B; 

%Error  by  every  Coordinate 
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Err  =  [abs  (Coorl  (1) -E  (1) ) , abs  (Coorl  (2) -E  (2) ) , 
abs  (Coorl  (3) -E  (3) ) ] ; 

%Error  Matrix 
Mat (i, : )  =  Err; 


end 


%calculate  Cov. Matrix 
Cov_Mat  =  cov (Mat); 

%Eigenvalues  and  eigenvectors 
[Eigenvectors ,  Eigenvalues  ]  =  eig  (Cov_Mat )  ; 

%Find  Diagonalized  Cov. Mat. 

Diagonal_C°v_Mat  = 

Eigen_Vectors . ' *inv (Cov_M&t ) *Eigen_Vectors ; 

%Std.  Dev  (includes  k  value) 

Std_Dev_Data  =  sqrt (diag (Eigenvalues) *k) ; 

Std_Dev_Data_Mat ( : , j )  =  Std_Dev_Data 

%Total  Standat  Deviation 
Tot_Std_Dev  = 

sqrt (Std_Dev_Data (1) A2+Std_Dev_Data (2) A2+Std_Dev_Data (3) A2) 


Tot_Std_Dev_Mat  (j,:)  =  Tot_Std_Dev 


end 

figure 

%Draw  the  error  ellipsoid  from  the  diagonalized  covariance 

matrix 

hold  on 

grid  on 

%equal  axis 
axis  equal 

%#  Change  the  camera  viewpoint 
view ([-36  18]); 

%label  axises  and  title 
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xlabel ( ' X-Axis ' ) 
ylabel ( ' Y-Axis ' ) 
zlabel ( ' Z-Axis ' ) 
title ( ' TDOA  Geometry') 

%plot  stationary  receivers 

plot3  (X2 , Y2 , Z2 ,  '  o  '  )  ; 
plot3  (X3 , Y3 , Z3 , ' o ' ) ; 
plot3  (X4 , Y4 , Z4 , 'o'  )  ; 
plot3  (X5 , Y5 , Z5 , ' o ' ) ; 

%plot  emitter 

plot3  (E  (1)  ,E  (2)  ,E  (3)  ,  'ms  '  )  ; 

%plot  moving  receiver 

plot 3 (XlCoor, YlCoor, ZlCoor, ' r — ' ) ; 

figure 

N  =  linspace ( 1 , m, m) ; 

%plot  the  result 

plot  (N, Tot_Std_Dev_Mat ) 

%label  axises  and  title 

xlabel (' Distance  Multiplier') 

ylabel  (' Total  Standard  Deviation  (m)  ') 

title  (' Effect  of  Distance  to  the  Accuracy  With  UAV') 
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MATLAB  CODE  FOR  SIMULATION  OF  CLOSED-FORM  GEOLOCATION 
TECHNIQUE  TO  SEE  THE  EFFECTS  OF  HAVING  TWO  MOVING  UAVS 
ALONG  WITH  STAIONARY  RECEIVERS  ON  ACCURACY. 


function  TDOA_Flying_two_UAV 
%  Written  by:  Volkan  TAS 

%  Based  on  code  developed  by  EZZAT  G.  BAKHOUM,  described  in 
%  Closed-Form  Solution  of  Hyperbolic  Geolocation 

%  Equations,  2006 

%  Latest  revision:  August  26,  2012 

%clear  the  screen 
clc 

%Variables 

%number  of  experiments  for  accuracy 
n  =  500; 

%number  of  the  positions  of  the  UAV 
m  =  2  0; 


k=8 . 934259824 ;  %k  for  95%  of  probability 


o 

o 

k=6. 922544695; 

%k 

for 

90% 

of 

probability 

o 

o 

k=5. 766298681; 

%k 

for 

85% 

of 

probability 

o, 

o 

k=4 . 949459443; 

%k 

for 

80% 

of 

probability 

o, 

o 

k=4 .314831079; 

%k 

for 

75% 

of 

probability 

%Coordinates  of  emitter 
E= [ 0  0  0  ]  ; 

%Speed  of  Light ( Propagation) 
c  =  3  *  1 0  A  8  ; 


%Set  All  Alphas  to  1 


alphal  =  1; 
alpha2  =  1; 
alpha3  =  1; 
alpha4  =  1; 
alpha5  =  1; 

%Stationary 

Receivers 

X3  =-1000; 

Y3  =  2000; 

Z3  = 

-2000; 

X4  =  -1000; 

Y4  =  -1000; 

Z4  = 

1000; 

X5  =  -3000; 

Y5  =  -1000; 

Z5  = 

-1000; 
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%for  loop  for  UAV  movement 
for  j  =  l:m, 

%locations  of  the  UAV 

XI  =  8  0  0  *  j ;  Y1  =  1 2  0  0  *  j ;  Z1  =  2000;  %Moving  UAV 

X2  =  7  0  0  *  j ;  Y2  =  -1800*j;  Z2  =  2000; 

%UAV  Coordinate  MAtrix 
XlCoor (1, j )  =  XI; 

YlCoor  (1, j )  =  Yl; 

ZlCoor ( 1 , j )  =  Zl; 

X2Coor  (1, j )  =  X2; 

Y2Coor (1, j )  =  Y2; 

Z2Coor ( 1 , j )  =  Z2 ; 

%for  loop  for  experiments  for  each  location  of  the  UAV 
for  i  =  1 : n, 

%Calculate  the  Distance  From  emitter  to  the 
Receiver  for  TOA 

ECl=sqrt (abs (  (E (1)  -XI) A2+ (E (2) -Yl) A2+  (E  (3) -Zl) A2) )  ; 
EC2=sqrt (abs (  (E (1) -X2) A2+ (E  (2) -Y2) A2+  (E  ( 3 ) -Z2 ) A2 ) )  ; 
EC3=sqrt (abs ( (E (1) -X3) A2+ (E (2) -Y3) A2+ (E (3) -Z3) A2) ) ; 
EC4=sqrt  (abs  (  (E  (1)  -X4)  A2+  (E  (2)  -Y4)  A2+(E(3)-Z4)A2)); 
EC5=sqrt  (abs  (  (E  (1)  -X5)  A2+  (E  (2)  -Y5)  A2+(E(3)-Z5)A2)); 

%Generate  random  noise  for  every  channel  related  to 
Standat  Deviation 

Noise  Err  =  random (' Normal 0 , 5e-8 , 1 , 5 ) ; 


%Calculate 

the  times 

of 

arrival 

. 

tl 

=  EC1 

/ 

(alphal 

* 

c) 

+ 

Noise 

Err 

(l)  ; 

t2 

=  EC2 

/ 

(alpha2 

* 

c) 

+ 

Noise 

Err (2) 

t3 

=  EC3 

/ 

(alpha3 

* 

c) 

+ 

Noise 

Err (3) ; 

1 4 

=  EC4 

/ 

(alpha4 

* 

c) 

+ 

Noise 

Err 

(4)  ; 

t5 

=  EC5 

/ 

(alpha5 

* 

c) 

+ 

Noise 

Err 

(5)  ; 

%from  now  on  calculate  the  emitter  location 

%Calculate  the  coefficients  in  Eqs . (17) .  Note  that 
alpha  (assumed) , not  alpha  actual  is  used. 

all  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2/(t3-tl)  * (X3/alpha3A2  -  XI /alphal A2 ) ; 
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al2  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2/(t3-tl)  * (Y3/alpha3A2  -  Y1 /alphal A2 ) ; 

al3  =  2/(t2-tl)  *  ( Z 2 / alpha2 A2  -  Zl/alphalA2)  - 

2/(t3-tl)  * (Z3/alpha3A2  -  Zl/alphalA2) ; 

bl  =  1/ ( t2-tl )  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t3-tl)  * 

( (X3A2+Y3A2+Z3A2) /alpha3A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( 1 3 - 1 2 ) ; 

a2 1  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (X4/alpha4A2  -  XI /alphal A2 ) ; 

a22  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2/(t4-tl)  * (Y4/alpha4A2  -  Y1 /alphal A2 ) ; 

a23  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2 /  ( 1 4-t 1 )  * (Z4/alpha4A2  -  Zl/alphalA2) ; 

b2  =  1/ (t2-tl)  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t4-tl)  * 

( (X4A2+Y4A2+Z4A2) /alpha4A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( 1 4 - 1 2 )  ; 

a31  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  - 
2/(t5-tl)  * (X5/alpha5A2  -  XI /alphal A2 ) ; 

a32  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  - 
2/(t5-tl)  * (Y5/alpha5A2  -  Y1 /alphal A2 ) ; 

a33  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  - 
2/(t5-tl)  * (Z5/alpha5A2  -  Zl/alphalA2) ; 

b3  =  1/ (t2-tl)  *  ( (X2A2+Y2A2+Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t5-tl)  * 

(  (X5A2+Y5A2  +  Z5A2) /alpha5A2-  (X1A2+Y1A2  +  Z1A2) /alphalA2)  + 
cA2  *  ( 1 5 - 1 2 )  ; 

%Calculate  the  matrix  A  on  the  left-hand  side. 

A  =  double ([[all  al2  al3] ;  [a21  a22  a23];  [a31  a32 

a  3  3  ]  ]  )  ; 


%Calculate  the  vector  B  on  the  right-hand  side. 
B  =  double ([bl;  b2;  b3]); 

%Solve  the  simultaneous  equations  AX  =  B. 

Coorl  =  (A' *A) \A ' *B; 

%Error  by  every  Coordinate 

Err  =  [abs (Coorl (1) -E  (1) ) , abs  (Coorl  (2) -E  (2) ) , 
abs  (Coorl  (3) -E  (3) ) ] ; 

%Error  Matrix 
Mat  (i,  : )  =  Err; 
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end 


%calculate  Cov. Matrix 
Cov_Mat  =  cov (Mat); 

%Eigenvalues  and  eigenvectors 
[Eigenvectors ,  Eigenvalues ]  =  eig  (Cov_Mat )  ; 

%Find  Diagonalized  Cov. Mat. 

Diagonal_C°v_^at  = 

Eigenvectors  _  1  *j_nv  (Cov_Mat )  *Eigen_Vectors ; 

%Std.  Dev  (includes  k  value) 

Std_Dev_Data  =  sqrt (diag (Eigenvalues) *k) ; 

Std_Dev_Data_Mat ( : , j )  =  Std_Dev_Data; 

%Total  Standat  Deviation 
Tot_Std_Dev  = 

sqrt (Std_Dev_Data (1) A2+Std_Dev_Data (2) A2+Std_Dev_Data (3) A2) 


Tot_Std_Dev_Mat  (j,:)  =  Tot_Std_Dev 


end 

figure 

%Draw  the  error  ellipsoid  from  the  diagonalized  covariance 

matrix 

hold  on 

grid  on 

%equal  axis 
axis  equal 

%#  Change  the  camera  viewpoint 
view ([-36  18]); 

%label  axises  and  title 
xlabel ( ' X-Axis ' ) 
ylabel ( ' Y-Axis ' ) 
zlabel ( ' Z-Axis ' ) 
title ( ' TDOA  Geometry') 
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%plot  stationary  receivers 


plot3  (X3, Y3, Z3, ' o ' ) ; 
plot3  (X4 , Y4 , Z4 , 'o'  )  ; 
plot3  (X5 , Y5 , Z5 , ' o ' ) ; 

%plot  emitter 

plot3  (E  (1)  ,E  (2)  ,E  (3)  ,  'ms  '  )  ; 

%plot  moving  receiver 

plot 3 (XlCoor, YlCoor, ZlCoor, ' r-- ' ) ; 

plot 3 (X2Coor , Y2Coor , Z2Coor , ' g — ' ) ; 

figure 

N  =  linspace ( 1 , m, m) ; 

%plot  the  result 

plot  (N, Tot_Std_Dev_Mat ) 

%label  axises  and  title 

xlabel (' Distance  Multiplier') 

ylabel  (' Total  Standard  Deviation  (m)  ') 

title  (' Effect  of  Distance  to  the  Accuracy  With  UAV') 
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o\°  o\° 


APPENDIX  G 


MATLAB  CODE  FOR  SIMULATION  OF  CLOSED-FORM  GEOLOCATION 
TECHNIQUE  TO  SEE  THE  EFFECTS  OF  THE  RECEIVERS  GEOMETRY  ON 

ACCURACY. 

function  TDOA_Receiver_Dist 
%  Written  by:  Volkan  TAS 

%  Based  on  code  developed  by  EZZAT  G.  BAKHOUM,  described  in 
%  Closed-Form  Solution  of  Hyperbolic  Geolocation 
%  Equations,  %  2006 
%  Latest  revision:  August  15,  2012 

%clear  the  screen 
clc 

%Variables 

%number  of  experiments  for  accuracy 
n  =  500; 

Mat  =  zeros (n, 3 ) ; 


k=8 . 934259824 ;  %k  for  95%  of  probability 


%  k=6. 922544695; 

%k 

for 

90% 

of 

probability 

%  k=5. 766298681; 

%k 

for 

85% 

of 

probability 

%  k=4 . 949459443; 

%k 

for 

80% 

of 

probability 

%  k=4 . 314831079; 

%k 

for 

75% 

of 

probability 

%Coordinates  of 

emitter 

E= [ 0  0  0  ]  ; 

%Speed  of  Light ( Propagation) 

c  =  3  *  1 0  A  8  ; 

%Set  All  Alphas  to  1 

alphal  =  1; 

alpha2  =  1; 

alpha3  =  1; 

alpha4  =  1; 

alpha5  =  1; 

Stationary  Receivers 
set  default  for  switch 
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val 


1; 


tmp  =  input  (sprintf  ('How  you  want  to  Distribute  the 
Receivers  \n  Around  the  Emitter  (Without  Altitude 
Difference) Press  1  \n  On  Staright  Line  (Without  Altitude 
Difference) Press  2  \n  Around  the  Emitter  (With  Altitude 
Difference)  Press  3  \n  On  Staright  Line  (With  Altitude 
Difference) Press  4  \n  Trapezoidal  Press  5  \n 
Parallelogram  Press  6  \n  Inverted  Triangle  Press  8 
' , 'val' ) ) ; 

%#  if  the  user  hits  'return'  without  writing  anything,  tmp 
is  empty  and  the  default  is  used 
if  -isempty ( tmp) 
val  =  tmp; 

end 

switch  val 


case  1 

%circle,  no 

elevation  difference 

XI 

=  2  e  3  ; 

Y1  = 

0; 

Zl  = 

1 . 9  7  e  3  ; 

X2 

=  7e2  ; 

Y2  = 

2e3; 

Z2  = 

1 . 9  8  e  3  ; 

X3 

=-2e3 ; 

Y3  = 

1 . 5  e  3 ; 

Z3  = 

1 . 96e3; 

X4 

=  -2e3; 

Y4  = 

- 1 . 5  e  3 ; 

Z4  = 

2e3; 

X5 

=  1  e  3  ; 

Y5  = 

-1 . 5e3 ; 

Z5  = 

1 . 99e3 ; 

case  2 

%0n  Staright 

Line 

(Without 

Altitude  Difference 

XI 

=  -2e3; 

Y1  = 

- 1  e  3  ; 

Zl  = 

0.99e3; 

X2 

=  - 1  e  3  ; 

Y2  = 

0; 

Z2  = 

le3  ; 

X3 

=  0; 

Y3  = 

2e3  ; 

Z3  = 

0 . 9  8  e  3  ; 

X4 

=  1  e  3  ; 

Y4  = 

3e3  ; 

Z4  = 

0. 97e3; 

X5 

=  1 . 5  e  3 ; 

Y5  = 

4e3 ; 

Z5  = 

1 . Ie3 ; 

case  3 

%Around  the 

Emitter  (With 

Altitude  Difference) 

XI 

=  2  e  3  ; 

Y1  = 

0; 

Zl  = 

0; 

X2 

=  0 . 7  e  3  ; 

Y2  = 

2e3  ; 

Z2  = 

le3  ; 

X3 

=-2e3 ; 

Y3  = 

1 . 5e3 ; 

Z3  = 

- 1  e  3  ; 

X4 

=  -2e3; 

Y4  = 

-1 . 5e3; 

Z4  = 

2e3  ; 

X5 

=  le3  ; 

Y5  = 

- 1 . 5  e  3 ; 

Z5  = 

-2e3  ; 

case  4 

%On  Staright 

Line 

(With  Altitude 

Difference) 

XI 

=  -2e3; 

Y1  = 

- 1  e  3  ; 

Zl  = 

0; 

X2 

=  - 1  e  3  ; 

Y2  = 

0; 

Z2  = 

le3  ; 

X3 

=  0; 

Y3  = 

2e3  ; 

Z3  = 

2e3  ; 

X4 

=  le3  ; 

Y4  = 

3e3  ; 

Z4  = 

3e3  ; 

X5 

=  1 . 5e3 ; 

Y5  = 

4e3  ; 

Z5  = 

4e3 ; 
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case  5  %Trapezoidal 


XI 

=  le3  ; 

Yl 

= 

0; 

Zl 

= 

0; 

X2 

=  2  e  3  ; 

Y2 

= 

2e3  ; 

Z2 

= 

-0.1e3; 

X3 

=  3e3; 

Y3 

= 

2e3  ; 

Z3 

= 

-0.1e3; 

X4 

=  4e3  ; 

Y4 

= 

2e3  ; 

Z4 

= 

-  0 . 1  e  3  ; 

X5 

=  5e3; 

Y5 

= 

0; 

Z5 

= 

-0.1e3; 

case  6 

%Parallelogram 

XI 

=  - 1  e  3  ; 

Yl 

= 

0; 

Zl 

= 

0; 

X2 

=  -2e3; 

Y2 

= 

2e3  ; 

Z2 

= 

-  0 . 1  e  3 ; 

X3 

=  0; 

Y3 

= 

2e3  ; 

Z3 

= 

-  0 . 1  e  3 ; 

X4 

=  -2e3; 

Y4 

= 

0; 

Z4 

= 

-  0 . 1  e  3 ; 

X5 

=  -3e3; 

Y5 

= 

0; 

Z5 

= 

-  0 . 1  e  3  ; 

case  7 

%Lorenze 

XI 

=  8e3; 

Yl 

= 

15e3; 

Zl 

= 

-  0 . 1 2  e  3  ; 

X2 

=  - 1 3  e  3  ; 

Y2 

= 

12.5e3; 

Z2 

= 

-0 . 13e3 ; 

X3 

=  0; 

Y3 

= 

1 8  e  3  ; 

Z3 

= 

-0.14e3; 

X4 

=  13e3 ; 

Y4 

= 

1 2 . 5  e  3 ; 

Z4 

= 

-  0 . 1 1  e  3 ; 

X5 

=  0; 

Y5 

= 

8; 

Z5 

= 

-  0 . 1  e  3 ; 

case  8 

%Inverted 

Triangle 

XI 

=  2  e  3  ; 

Yl 

= 

0; 

Zl 

= 

0; 

X2 

=  4  e  3  ; 

Y2 

= 

2e3  ; 

Z2 

= 

-  0 . 1  e  3 ; 

X3 

=  6e3; 

Y3 

= 

-2e3; 

Z3 

= 

-0.1e3; 

X4 

=  8e3; 

Y4 

= 

4e3; 

Z4 

= 

-0.1e3; 

X5 

=  1 0  e  3 ; 

Y5 

= 

-  4  e  3 ; 

Z5 

= 

-  0 . 1  e  3 ; 

end 

for  i  =  1 : n, 

%Calculate  the  Distance  From  emitter  to  the  Receiver 
for  TOA 

ECl=sqrt (abs (  (E (1) -XI) A2+ (E (2) -Yl) A2  + (E  (3) -Zl) A2) )  ; 
EC2=sqrt (abs (  (E (1) -X2) A2+ (E (2) -Y2) A2+  (E  ( 3 ) -Z2 ) A2 ) )  ; 
EC3=sqrt (abs ( (E (1) -X3) A2+ (E (2) -Y3) A2+ (E (3) -Z3) A2) ) ; 
EC4=sqrt (abs (  (E (1) -X4) A2+  (E  (2) -Y4) A2+  (E  ( 3 ) -Z4 ) A2 ) )  ; 
EC5=sqrt (abs (  (E (1) -X5) A2+ (E (2) -Y5) A2+  (E  (3) -Z5) A2) )  ; 

%Generate  random  noise  for  every  channel  related  to 
Standat  Deviation 

Noise_Err  =  random (' Normal 0 , 5e-8 , 1 , 5 ) ; 

%Calculate  the  times  of  arrival. 

tl  =  EC1  /  (alphal  *  c)  +  Noise_Err ( 1 ) ; 
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t2 

=  EC2 

/ 

(alpha2 

* 

c) 

+ 

Noise 

Err  (2)  ; 

t3 

=  EC3 

/ 

(alpha3 

* 

c) 

+ 

Noise 

_Err  (3)  ; 

1 4 

=  EC4 

/ 

(alpha4 

* 

c) 

+ 

Noise 

Err ( 4 ) ; 

t5 

=  EC5 

/ 

(alpha5 

* 

c) 

+ 

Noise 

Err ( 5 ) ; 

%from  now  on  calculate  the  emitter  location 


%Calculate  the  coefficients  in  Eqs . (17) .  Note  that 
alpha  (assumed) , not  alpha  actual  is  used. 

all  =  2/(t2-tl)  *  (X2/alpha2A2  -  Xl/alphalA2)  -  2/ (t3 
tl)  * (X3/alpha3A2  -  XI /alphal A2 ) ; 

al2  =  2 /  ( t2 - 1 1 )  *  (Y2/alpha2A2  -  Yl/alphalA2)  -  2/ (t3 
tl)  * (Y3/alpha3A2  -  Y1 /alphal A2 ) ; 

al3  =  2 /  ( t2 - 1 1 )  *  (Z2/alpha2A2  -  Zl/alphalA2)  -  2/ (t3 
tl)  * (Z3/alpha3A2  -  Zl/alphalA2) ; 

bl  =  1/ (t2-tl)  *  (  (X2A2+Y2A2  +  Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t3-tl)  * 

(  (X3A2+Y3A2  +  Z3A2) /alpha3A2-  (X1A2+Y1A2  +  Z1A2) /alphalA2)  + 
cA2  *  ( t3-t2 )  ; 

a2 1  =  2 /  ( t2 - 1 1 )  *  (X2/alpha2A2  -  Xl/alphalA2)  -  2/ (t4 
tl)  * (X4/alpha4A2  -  XI /alphal A2 ) ; 

a22  =  2/(t2-tl)  *  (Y2/alpha2A2  -  Yl/alphalA2)  -  2/ (t4 
tl)  * ( Y4 /alpha4 A2  -  Y1 /alphal A2 ) ; 

a23  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  -  2/ (t4 
tl)  * (Z4/alpha4A2  -  Zl/alphalA2) ; 

b2  =  1 /  ( t2-t 1 )  *  ( (X2A2+Y2A2  +  Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t4-tl)  * 

( (X4A2+Y4A2+Z4A2) /alpha4A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( 1 4 - 1 2 )  ; 

a31  =  2 /  ( t2 - 1 1 )  *  (X2/alpha2A2  -  Xl/alphalA2)  -  2/ (t5 
tl)  * (X5/alpha5A2  -  XI /alphal A2 ) ; 

a32  =  2 /  ( t2 - 1 1 )  *  (Y2/alpha2A2  -  Yl/alphalA2)  -  2/ (t5 
tl)  * (Y5/alpha5A2  -  Y1 /alphal A2 ) ; 

a33  =  2/(t2-tl)  *  (Z2/alpha2A2  -  Zl/alphalA2)  -  2/ (t5 
tl)  * (Z5/alpha5A2  -  Zl/alphalA2) ; 

b3  =  1 /  ( t2-t 1 )  *  ( (X2A2+Y2A2  +  Z2A2) /alpha2A2- 

(XI A2+Y1 A2+Z 1 A2 ) / alphal A 2 ) -  l/(t5-tl)  * 

( (X5A2+Y5A2+Z5A2) /alpha5A2-  (X1A2+Y1A2+Z1A2) /alphalA2)  + 
cA2  *  ( t5-t2 )  ; 

%Calculate  the  matrix  A  on  the  left-hand  side. 

A  =  double ([[all  al2  al3] ;  [a21  a22  a23];  [a31  a32 

a  3  3  ]  ]  )  ; 
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%Calculate  the  vector  B  on  the  right-hand  side. 
B  =  double ( [bl;  b2;  b3]); 

%Solve  the  simultaneous  equations  AX  =  B. 

Coor  =  (A' *A) \A ' *B; 


%Error  by  every  Coordinate 

Err  =  [abs (Coor  (1) -E (1) ) , abs  (Coor  (2) -E  (2) ) , 
abs  (Coor  (3) -E  (3) ) ] ; 

%Error  Matrix 
Mat  (i,  : )  =  Err; 


end 

%calculate  Cov. Matrix 
Cov_Mat  =  cov (Mat); 

%Eigenvalues  and  eigenvectors 
[Eigenvectors ,  Eigenvalues  ]  =  eig (Cov_Mat ) ; 

%Find  angles  (from  "Computing  Euler  angles  from  a  rotation 
matrix" ) 

Beta  =  -asind  (Eigenvectors  ( 3 ,  1 )  )  ; 

Gamma  = 

(atan2  (Eigen_Vectors (2, 1) /cosd(Beta) , Eigen_Vectors ( 1 ,  1 )  /  cos 
d(Beta) ))* (180/pi) ; 

Alpha  = 

(atan2 (Eigen_Vectors (3,2) /cosd(Beta) , Eigen_Vectors (3, 3) /cos 
d(Beta) ))* (180/pi) ; 

%Std.  Dev  (includes  k  value) 

Std_Dev_Data  =  sqrt (diag (Eigenvalues) *k) 

%Total  Standat  Deviation 
Tot_Std_Dev  = 

sqrt (Std_Dev_Data (1) A2+Std_Dev_Data (2) A2+Std_Dev_Data (3) A2) 

%Draw  the  error  ellipsoid  from  the  diagonalized  covariance 

matrix 

hold  on 

grid  on 
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%ellipsoid 
[x,  y,  z]  = 

ellipsoid(E (1) , E (2) , E (3) , Std_Dev_Data (1) , Std_Dev_Data (2) , St 
d_Dev_Data ( 3 ) ,20) ; 

%color  of  the  ellisoid 
colormap  copper; 

%ploy  ellipsoid 
hMesh  =  mesh (x, y, z) ; 

%rotate  ellipsoid 
rotate (hMesh, [ 1  0  0], Alpha); 
rotate (hMesh, [ 0  1  0],Beta); 
rotate (hMesh, [0  0  1 ], Gamma); 

%equal  axis 
axis  equal 

%#  Change  the  camera  viewpoint 
view ([-36  18]); 

%label  axises  and  title 
xlabel ( ' X-Axis ' ) 
ylabel ( ' Y-Axis ' ) 
zlabel ( ' Z-Axis ' ) 

title ( ' TDOA  Geometry  and  Confidence  Ellipsoid') 

%plot  stationary  receivers 
plot3  (XI, Yl, Zl, 'o'  )  ; 
plot3  (X2 , Y2 , Z2 , 'o' ) ; 
plot3  (X3, Y3, Z3, 'o' ) ; 
plot3  (X4 , Y4 , Z4 , 'o' ) ; 
plot3  (X5 , Y5 , Z5 , ' o ' ) ; 

%plot  emitter 

plot3  (E  (1) ,E  (2) ,E  (3) ,  'ms ' ) ; 
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